home *** CD-ROM | disk | FTP | other *** search
/ NeXTSTEP 3.3 (Developer)…68k, x86, SPARC, PA-RISC] / NeXTSTEP 3.3 Dev Intel.iso / NextDeveloper / Source / GNU / cc / real.c < prev    next >
C/C++ Source or Header  |  1993-10-13  |  118KB  |  6,053 lines
  1. /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
  2. and support for XFmode IEEE extended real floating point arithmetic.
  3. Contributed by Stephen L. Moshier (moshier@world.std.com).
  4.  
  5.    Copyright (C) 1993 Free Software Foundation, Inc.
  6.  
  7. This file is part of GNU CC.
  8.  
  9. GNU CC is free software; you can redistribute it and/or modify
  10. it under the terms of the GNU General Public License as published by
  11. the Free Software Foundation; either version 2, or (at your option)
  12. any later version.
  13.  
  14. GNU CC is distributed in the hope that it will be useful,
  15. but WITHOUT ANY WARRANTY; without even the implied warranty of
  16. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17. GNU General Public License for more details.
  18.  
  19. You should have received a copy of the GNU General Public License
  20. along with GNU CC; see the file COPYING.  If not, write to
  21. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  22.  
  23. #include <stdio.h>
  24. #include <errno.h>
  25. #include "config.h"
  26. #include "tree.h"
  27.  
  28. #ifndef errno
  29. extern int errno;
  30. #endif
  31.  
  32. /* To enable support of XFmode extended real floating point, define
  33. LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  34.  
  35. To support cross compilation between IEEE, VAX and IBM floating
  36. point formats, define REAL_ARITHMETIC in the tm.h file.
  37.  
  38. In either case the machine files (tm.h) must not contain any code
  39. that tries to use host floating point arithmetic to convert
  40. REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  41. etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  42. be intelligible to the host computer's native arithmetic.
  43.  
  44. The emulator defaults to the host's floating point format so that
  45. its decimal conversion functions can be used if desired (see
  46. real.h).
  47.  
  48. The first part of this file interfaces gcc to ieee.c, which is a
  49. floating point arithmetic suite that was not written with gcc in
  50. mind.  The interface is followed by ieee.c itself and related
  51. items. Avoid changing ieee.c unless you have suitable test
  52. programs available.  A special version of the PARANOIA floating
  53. point arithmetic tester, modified for this purpose, can be found
  54. on usc.edu : /pub/C-numanal/ieeetest.zoo.  Some tutorial
  55. information on ieee.c is given in my book: S. L. Moshier,
  56. _Methods and Programs for Mathematical Functions_, Prentice-Hall
  57. or Simon & Schuster Int'l, 1989.  A library of XFmode elementary
  58. transcendental functions can be obtained by ftp from
  59. research.att.com: netlib/cephes/ldouble.shar.Z  */
  60.  
  61. /* Type of computer arithmetic.
  62.  * Only one of DEC, IBM, MIEEE, IBMPC, or UNK should get defined.
  63.  */
  64.  
  65. /* `MIEEE' refers generically to big-endian IEEE floating-point data
  66.    structure.  This definition should work in SFmode `float' type and
  67.    DFmode `double' type on virtually all big-endian IEEE machines.
  68.    If LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then MIEEE
  69.    also invokes the particular XFmode (`long double' type) data
  70.    structure used by the Motorola 680x0 series processors.
  71.  
  72.    `IBMPC' refers generally to little-endian IEEE machines. In this
  73.    case, if LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then
  74.    IBMPC also invokes the particular XFmode `long double' data
  75.    structure used by the Intel 80x86 series processors.
  76.  
  77.    `DEC' refers specifically to the Digital Equipment Corp PDP-11
  78.    and VAX floating point data structure.  This model currently
  79.    supports no type wider than DFmode.
  80.  
  81.    `IBM' refers specifically to the IBM System/370 and compatible
  82.    floating point data structure.  This model currently supports
  83.    no type wider than DFmode.  The IBM conversions were contributed by
  84.    frank@atom.ansto.gov.au (Frank Crawford).
  85.  
  86.    If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
  87.    then `long double' and `double' are both implemented, but they
  88.    both mean DFmode.  In this case, the software floating-point
  89.    support available here is activated by writing
  90.       #define REAL_ARITHMETIC
  91.    in tm.h. 
  92.  
  93.    The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
  94.    and may deactivate XFmode since `long double' is used to refer
  95.    to both modes.
  96.  
  97.    The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
  98.    contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
  99.    separate the floating point unit's endian-ness from that of
  100.    the integer addressing.  This permits one to define a big-endian
  101.    FPU on a little-endian machine (e.g., ARM).  An extension to
  102.    BYTES_BIG_ENDIAN may be required for some machines in the future.
  103.    These optional macros may be defined in tm.h.  In real.h, they
  104.    default to WORDS_BIG_ENDIAN, etc., so there is no need to define
  105.    them for any normal host or target machine on which the floats
  106.    and the integers have the same endian-ness.   */
  107.  
  108.  
  109. /* The following converts gcc macros into the ones used by this file.  */
  110.  
  111. /* REAL_ARITHMETIC defined means that macros in real.h are
  112.    defined to call emulator functions.  */
  113. #ifdef REAL_ARITHMETIC
  114.  
  115. #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  116. /* PDP-11, Pro350, VAX: */
  117. #define DEC 1
  118. #else /* it's not VAX */
  119. #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  120. /* IBM System/370 style */
  121. #define IBM 1
  122. #else /* it's also not an IBM */
  123. #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  124. #if FLOAT_WORDS_BIG_ENDIAN
  125. /* Motorola IEEE, high order words come first (Sun workstation): */
  126. #define MIEEE 1
  127. #else /* not big-endian */
  128. /* Intel IEEE, low order words come first:
  129.  */
  130. #define IBMPC 1
  131. #endif /*  big-endian */
  132. #else /* it's not IEEE either */
  133. /* UNKnown arithmetic.  We don't support this and can't go on. */
  134. unknown arithmetic type
  135. #define UNK 1
  136. #endif /* not IEEE */
  137. #endif /* not IBM */
  138. #endif /* not VAX */
  139.  
  140. #else
  141. /* REAL_ARITHMETIC not defined means that the *host's* data
  142.    structure will be used.  It may differ by endian-ness from the
  143.    target machine's structure and will get its ends swapped
  144.    accordingly (but not here).  Probably only the decimal <-> binary
  145.    functions in this file will actually be used in this case.  */
  146. #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  147. #define DEC 1
  148. #else /* it's not VAX */
  149. #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  150. /* IBM System/370 style */
  151. #define IBM 1
  152. #else /* it's also not an IBM */
  153. #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  154. #if HOST_FLOAT_WORDS_BIG_ENDIAN
  155. #define MIEEE 1
  156. #else /* not big-endian */
  157. #define IBMPC 1
  158. #endif /*  big-endian */
  159. #else /* it's not IEEE either */
  160. unknown arithmetic type
  161. #define UNK 1
  162. #endif /* not IEEE */
  163. #endif /* not IBM */
  164. #endif /* not VAX */
  165.  
  166. #endif /* REAL_ARITHMETIC not defined */
  167.  
  168. /* Define INFINITY for support of infinity.
  169.    Define NANS for support of Not-a-Number's (NaN's).  */
  170. #if !defined(DEC) && !defined(IBM)
  171. #define INFINITY
  172. #define NANS
  173. #endif
  174.  
  175. /* Support of NaNs requires support of infinity. */
  176. #ifdef NANS
  177. #ifndef INFINITY
  178. #define INFINITY
  179. #endif
  180. #endif
  181.  
  182. /* Find a host integer type that is at least 16 bits wide,
  183.    and another type at least twice whatever that size is. */
  184.  
  185. #if HOST_BITS_PER_CHAR >= 16
  186. #define EMUSHORT char
  187. #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
  188. #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
  189. #else
  190. #if HOST_BITS_PER_SHORT >= 16
  191. #define EMUSHORT short
  192. #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
  193. #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
  194. #else
  195. #if HOST_BITS_PER_INT >= 16
  196. #define EMUSHORT int
  197. #define EMUSHORT_SIZE HOST_BITS_PER_INT
  198. #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
  199. #else
  200. #if HOST_BITS_PER_LONG >= 16
  201. #define EMUSHORT long
  202. #define EMUSHORT_SIZE HOST_BITS_PER_LONG
  203. #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
  204. #else
  205. /*  You will have to modify this program to have a smaller unit size. */
  206. #define EMU_NON_COMPILE
  207. #endif
  208. #endif
  209. #endif
  210. #endif
  211.  
  212. #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
  213. #define EMULONG short
  214. #else
  215. #if HOST_BITS_PER_INT >= EMULONG_SIZE
  216. #define EMULONG int
  217. #else
  218. #if HOST_BITS_PER_LONG >= EMULONG_SIZE
  219. #define EMULONG long
  220. #else
  221. #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
  222. #define EMULONG long long int
  223. #else
  224. /*  You will have to modify this program to have a smaller unit size. */
  225. #define EMU_NON_COMPILE
  226. #endif
  227. #endif
  228. #endif
  229. #endif
  230.  
  231.  
  232. /* The host interface doesn't work if no 16-bit size exists. */
  233. #if EMUSHORT_SIZE != 16
  234. #define EMU_NON_COMPILE
  235. #endif
  236.  
  237. /* OK to continue compilation. */
  238. #ifndef EMU_NON_COMPILE
  239.  
  240. /* Construct macros to translate between REAL_VALUE_TYPE and e type.
  241.    In GET_REAL and PUT_REAL, r and e are pointers.
  242.    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
  243.    in memory, with no holes.  */
  244.  
  245. #if LONG_DOUBLE_TYPE_SIZE == 96
  246. /* Number of 16 bit words in external e type format */
  247. #define NE 6
  248. #define MAXDECEXP 4932
  249. #define MINDECEXP -4956
  250. #define GET_REAL(r,e) bcopy (r, e, 2*NE)
  251. #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
  252. #else /* no XFmode */
  253. #if LONG_DOUBLE_TYPE_SIZE == 128
  254. #define NE 10
  255. #define MAXDECEXP 4932
  256. #define MINDECEXP -4977
  257. #define GET_REAL(r,e) bcopy (r, e, 2*NE)
  258. #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
  259. #else
  260. #define NE 6
  261. #define MAXDECEXP 4932
  262. #define MINDECEXP -4956
  263. #ifdef REAL_ARITHMETIC
  264. /* Emulator uses target format internally
  265.    but host stores it in host endian-ness. */
  266.  
  267. #if HOST_FLOAT_WORDS_BIG_ENDIAN == FLOAT_WORDS_BIG_ENDIAN
  268. #define GET_REAL(r,e) e53toe ((r), (e))
  269. #define PUT_REAL(e,r) etoe53 ((e), (r))
  270.  
  271. #else /* endian-ness differs */
  272. /* emulator uses target endian-ness internally */
  273. #define GET_REAL(r,e)        \
  274. do { EMUSHORT w[4];        \
  275.  w[3] = ((EMUSHORT *) r)[0];    \
  276.  w[2] = ((EMUSHORT *) r)[1];    \
  277.  w[1] = ((EMUSHORT *) r)[2];    \
  278.  w[0] = ((EMUSHORT *) r)[3];    \
  279.  e53toe (w, (e)); } while (0)
  280.  
  281. #define PUT_REAL(e,r)        \
  282. do { EMUSHORT w[4];        \
  283.  etoe53 ((e), w);        \
  284.  *((EMUSHORT *) r) = w[3];    \
  285.  *((EMUSHORT *) r + 1) = w[2];    \
  286.  *((EMUSHORT *) r + 2) = w[1];    \
  287.  *((EMUSHORT *) r + 3) = w[0]; } while (0)
  288.  
  289. #endif /* endian-ness differs */
  290.  
  291. #else /* not REAL_ARITHMETIC */
  292.  
  293. /* emulator uses host format */
  294. #define GET_REAL(r,e) e53toe ((r), (e))
  295. #define PUT_REAL(e,r) etoe53 ((e), (r))
  296.  
  297. #endif /* not REAL_ARITHMETIC */
  298. #endif /* not TFmode */
  299. #endif /* no XFmode */
  300.  
  301.  
  302. /* Number of 16 bit words in internal format */
  303. #define NI (NE+3)
  304.  
  305. /* Array offset to exponent */
  306. #define E 1
  307.  
  308. /* Array offset to high guard word */
  309. #define M 2
  310.  
  311. /* Number of bits of precision */
  312. #define NBITS ((NI-4)*16)
  313.  
  314. /* Maximum number of decimal digits in ASCII conversion
  315.  * = NBITS*log10(2)
  316.  */
  317. #define NDEC (NBITS*8/27)
  318.  
  319. /* The exponent of 1.0 */
  320. #define EXONE (0x3fff)
  321.  
  322. void warning ();
  323. extern int extra_warnings;
  324. int ecmp (), enormlz (), eshift ();
  325. int eisneg (), eisinf (), eisnan (), eiisinf (), eiisnan ();
  326. void eadd (), esub (), emul (), ediv ();
  327. void eshup1 (), eshup8 (), eshup6 (), eshdn1 (), eshdn8 (), eshdn6 ();
  328. void eabs (), eneg (), emov (), eclear (), einfin (), efloor ();
  329. void eldexp (), efrexp (), eifrac (), euifrac (), ltoe (), ultoe ();
  330. void ereal_to_decimal (), eiinfin (), einan ();
  331. void esqrt (), elog (), eexp (), etanh (), epow ();
  332. void asctoe (), asctoe24 (), asctoe53 (), asctoe64 (), asctoe113 ();
  333. void etoasc (), e24toasc (), e53toasc (), e64toasc (), e113toasc ();
  334. void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
  335. void etoe113 (), e113toe ();
  336. void mtherr (), make_nan ();
  337. void enan ();
  338. extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
  339. extern unsigned EMUSHORT elog2[], esqrt2[];
  340.  
  341. /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
  342.    swapping ends if required, into output array of longs.  The
  343.    result is normally passed to fprintf by the ASM_OUTPUT_ macros.   */
  344. void 
  345. endian (e, x, mode)
  346.      unsigned EMUSHORT e[];
  347.      long x[];
  348.      enum machine_mode mode;
  349. {
  350.   unsigned long th, t;
  351.  
  352. #if FLOAT_WORDS_BIG_ENDIAN
  353.   switch (mode)
  354.     {
  355.  
  356.     case TFmode:
  357.       /* Swap halfwords in the fourth long. */
  358.       th = (unsigned long) e[6] & 0xffff;
  359.       t = (unsigned long) e[7] & 0xffff;
  360.       t |= th << 16;
  361.       x[3] = (long) t;
  362.  
  363.     case XFmode:
  364.  
  365.       /* Swap halfwords in the third long. */
  366.       th = (unsigned long) e[4] & 0xffff;
  367.       t = (unsigned long) e[5] & 0xffff;
  368.       t |= th << 16;
  369.       x[2] = (long) t;
  370.       /* fall into the double case */
  371.  
  372.     case DFmode:
  373.  
  374.       /* swap halfwords in the second word */
  375.       th = (unsigned long) e[2] & 0xffff;
  376.       t = (unsigned long) e[3] & 0xffff;
  377.       t |= th << 16;
  378.       x[1] = (long) t;
  379.       /* fall into the float case */
  380.  
  381.     case SFmode:
  382.  
  383.       /* swap halfwords in the first word */
  384.       th = (unsigned long) e[0] & 0xffff;
  385.       t = (unsigned long) e[1] & 0xffff;
  386.       t |= th << 16;
  387.       x[0] = t;
  388.       break;
  389.  
  390.     default:
  391.       abort ();
  392.     }
  393.  
  394. #else
  395.  
  396.   /* Pack the output array without swapping. */
  397.  
  398.   switch (mode)
  399.     {
  400.  
  401.     case TFmode:
  402.  
  403.       /* Pack the fourth long. */
  404.       th = (unsigned long) e[7] & 0xffff;
  405.       t = (unsigned long) e[6] & 0xffff;
  406.       t |= th << 16;
  407.       x[3] = (long) t;
  408.  
  409.     case XFmode:
  410.  
  411.       /* Pack the third long.
  412.      Each element of the input REAL_VALUE_TYPE array has 16 useful bits
  413.      in it.  */
  414.       th = (unsigned long) e[5] & 0xffff;
  415.       t = (unsigned long) e[4] & 0xffff;
  416.       t |= th << 16;
  417.       x[2] = (long) t;
  418.       /* fall into the double case */
  419.  
  420.     case DFmode:
  421.  
  422.       /* pack the second long */
  423.       th = (unsigned long) e[3] & 0xffff;
  424.       t = (unsigned long) e[2] & 0xffff;
  425.       t |= th << 16;
  426.       x[1] = (long) t;
  427.       /* fall into the float case */
  428.  
  429.     case SFmode:
  430.  
  431.       /* pack the first long */
  432.       th = (unsigned long) e[1] & 0xffff;
  433.       t = (unsigned long) e[0] & 0xffff;
  434.       t |= th << 16;
  435.       x[0] = t;
  436.       break;
  437.  
  438.     default:
  439.       abort ();
  440.     }
  441.  
  442. #endif
  443. }
  444.  
  445.  
  446. /* This is the implementation of the REAL_ARITHMETIC macro.
  447.  */
  448. void 
  449. earith (value, icode, r1, r2)
  450.      REAL_VALUE_TYPE *value;
  451.      int icode;
  452.      REAL_VALUE_TYPE *r1;
  453.      REAL_VALUE_TYPE *r2;
  454. {
  455.   unsigned EMUSHORT d1[NE], d2[NE], v[NE];
  456.   enum tree_code code;
  457.  
  458.   GET_REAL (r1, d1);
  459.   GET_REAL (r2, d2);
  460. #ifdef NANS
  461. /*  Return NaN input back to the caller. */
  462.   if (eisnan (d1))
  463.     {
  464.       PUT_REAL (d1, value);
  465.       return;
  466.     }
  467.   if (eisnan (d2))
  468.     {
  469.       PUT_REAL (d2, value);
  470.       return;
  471.     }
  472. #endif
  473.   code = (enum tree_code) icode;
  474.   switch (code)
  475.     {
  476.     case PLUS_EXPR:
  477.       eadd (d2, d1, v);
  478.       break;
  479.  
  480.     case MINUS_EXPR:
  481.       esub (d2, d1, v);        /* d1 - d2 */
  482.       break;
  483.  
  484.     case MULT_EXPR:
  485.       emul (d2, d1, v);
  486.       break;
  487.  
  488.     case RDIV_EXPR:
  489. #ifndef REAL_INFINITY
  490.       if (ecmp (d2, ezero) == 0)
  491.     {
  492. #ifdef NANS
  493.     enan (v);
  494.     break;
  495. #else
  496.     abort ();
  497. #endif
  498.     }
  499. #endif
  500.       ediv (d2, d1, v);    /* d1/d2 */
  501.       break;
  502.  
  503.     case MIN_EXPR:        /* min (d1,d2) */
  504.       if (ecmp (d1, d2) < 0)
  505.     emov (d1, v);
  506.       else
  507.     emov (d2, v);
  508.       break;
  509.  
  510.     case MAX_EXPR:        /* max (d1,d2) */
  511.       if (ecmp (d1, d2) > 0)
  512.     emov (d1, v);
  513.       else
  514.     emov (d2, v);
  515.       break;
  516.     default:
  517.       emov (ezero, v);
  518.       break;
  519.     }
  520. PUT_REAL (v, value);
  521. }
  522.  
  523.  
  524. /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
  525.  * implements REAL_VALUE_RNDZINT (x) (etrunci (x))
  526.  */
  527. REAL_VALUE_TYPE 
  528. etrunci (x)
  529.      REAL_VALUE_TYPE x;
  530. {
  531.   unsigned EMUSHORT f[NE], g[NE];
  532.   REAL_VALUE_TYPE r;
  533.   HOST_WIDE_INT l;
  534.  
  535.   GET_REAL (&x, g);
  536. #ifdef NANS
  537.   if (eisnan (g))
  538.     return (x);
  539. #endif
  540.   eifrac (g, &l, f);
  541.   ltoe (&l, g);
  542.   PUT_REAL (g, &r);
  543.   return (r);
  544. }
  545.  
  546.  
  547. /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
  548.  * implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x))
  549.  */
  550. REAL_VALUE_TYPE 
  551. etruncui (x)
  552.      REAL_VALUE_TYPE x;
  553. {
  554.   unsigned EMUSHORT f[NE], g[NE];
  555.   REAL_VALUE_TYPE r;
  556.   unsigned HOST_WIDE_INT l;
  557.  
  558.   GET_REAL (&x, g);
  559. #ifdef NANS
  560.   if (eisnan (g))
  561.     return (x);
  562. #endif
  563.   euifrac (g, &l, f);
  564.   ultoe (&l, g);
  565.   PUT_REAL (g, &r);
  566.   return (r);
  567. }
  568.  
  569.  
  570. /* This is the REAL_VALUE_ATOF function.
  571.  * It converts a decimal string to binary, rounding off
  572.  * as indicated by the machine_mode argument.  Then it
  573.  * promotes the rounded value to REAL_VALUE_TYPE.
  574.  */
  575. REAL_VALUE_TYPE 
  576. ereal_atof (s, t)
  577.      char *s;
  578.      enum machine_mode t;
  579. {
  580.   unsigned EMUSHORT tem[NE], e[NE];
  581.   REAL_VALUE_TYPE r;
  582.  
  583.   switch (t)
  584.     {
  585.     case SFmode:
  586.       asctoe24 (s, tem);
  587.       e24toe (tem, e);
  588.       break;
  589.     case DFmode:
  590.       asctoe53 (s, tem);
  591.       e53toe (tem, e);
  592.       break;
  593.     case XFmode:
  594.       asctoe64 (s, tem);
  595.       e64toe (tem, e);
  596.       break;
  597.     case TFmode:
  598.       asctoe113 (s, tem);
  599.       e113toe (tem, e);
  600.       break;
  601.     default:
  602.       asctoe (s, e);
  603.     }
  604.   PUT_REAL (e, &r);
  605.   return (r);
  606. }
  607.  
  608.  
  609. /* Expansion of REAL_NEGATE.
  610.  */
  611. REAL_VALUE_TYPE 
  612. ereal_negate (x)
  613.      REAL_VALUE_TYPE x;
  614. {
  615.   unsigned EMUSHORT e[NE];
  616.   REAL_VALUE_TYPE r;
  617.  
  618.   GET_REAL (&x, e);
  619. #ifdef NANS
  620.   if (eisnan (e))
  621.     return (x);
  622. #endif
  623.   eneg (e);
  624.   PUT_REAL (e, &r);
  625.   return (r);
  626. }
  627.  
  628.  
  629. /* Round real toward zero to HOST_WIDE_INT
  630.  * implements REAL_VALUE_FIX (x).
  631.  */
  632. HOST_WIDE_INT
  633. efixi (x)
  634.      REAL_VALUE_TYPE x;
  635. {
  636.   unsigned EMUSHORT f[NE], g[NE];
  637.   HOST_WIDE_INT l;
  638.  
  639.   GET_REAL (&x, f);
  640. #ifdef NANS
  641.   if (eisnan (f))
  642.     {
  643.       warning ("conversion from NaN to int");
  644.       return (-1);
  645.     }
  646. #endif
  647.   eifrac (f, &l, g);
  648.   return l;
  649. }
  650.  
  651. /* Round real toward zero to unsigned HOST_WIDE_INT
  652.  * implements  REAL_VALUE_UNSIGNED_FIX (x).
  653.  * Negative input returns zero.
  654.  */
  655. unsigned HOST_WIDE_INT
  656. efixui (x)
  657.      REAL_VALUE_TYPE x;
  658. {
  659.   unsigned EMUSHORT f[NE], g[NE];
  660.   unsigned HOST_WIDE_INT l;
  661.  
  662.   GET_REAL (&x, f);
  663. #ifdef NANS
  664.   if (eisnan (f))
  665.     {
  666.       warning ("conversion from NaN to unsigned int");
  667.       return (-1);
  668.     }
  669. #endif
  670.   euifrac (f, &l, g);
  671.   return l;
  672. }
  673.  
  674.  
  675. /* REAL_VALUE_FROM_INT macro.
  676.  */
  677. void 
  678. ereal_from_int (d, i, j)
  679.      REAL_VALUE_TYPE *d;
  680.      HOST_WIDE_INT i, j;
  681. {
  682.   unsigned EMUSHORT df[NE], dg[NE];
  683.   HOST_WIDE_INT low, high;
  684.   int sign;
  685.  
  686.   sign = 0;
  687.   low = i;
  688.   if ((high = j) < 0)
  689.     {
  690.       sign = 1;
  691.       /* complement and add 1 */
  692.       high = ~high;
  693.       if (low)
  694.     low = -low;
  695.       else
  696.     high += 1;
  697.     }
  698.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  699.   ultoe (&high, dg);
  700.   emul (dg, df, dg);
  701.   ultoe (&low, df);
  702.   eadd (df, dg, dg);
  703.   if (sign)
  704.     eneg (dg);
  705.   PUT_REAL (dg, d);
  706. }
  707.  
  708.  
  709. /* REAL_VALUE_FROM_UNSIGNED_INT macro.
  710.  */
  711. void 
  712. ereal_from_uint (d, i, j)
  713.      REAL_VALUE_TYPE *d;
  714.      unsigned HOST_WIDE_INT i, j;
  715. {
  716.   unsigned EMUSHORT df[NE], dg[NE];
  717.   unsigned HOST_WIDE_INT low, high;
  718.  
  719.   low = i;
  720.   high = j;
  721.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  722.   ultoe (&high, dg);
  723.   emul (dg, df, dg);
  724.   ultoe (&low, df);
  725.   eadd (df, dg, dg);
  726.   PUT_REAL (dg, d);
  727. }
  728.  
  729.  
  730. /* REAL_VALUE_TO_INT macro
  731.  */
  732. void 
  733. ereal_to_int (low, high, rr)
  734.      HOST_WIDE_INT *low, *high;
  735.      REAL_VALUE_TYPE rr;
  736. {
  737.   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
  738.   int s;
  739.  
  740.   GET_REAL (&rr, d);
  741. #ifdef NANS
  742.   if (eisnan (d))
  743.     {
  744.       warning ("conversion from NaN to int");
  745.       *low = -1;
  746.       *high = -1;
  747.       return;
  748.     }
  749. #endif
  750.   /* convert positive value */
  751.   s = 0;
  752.   if (eisneg (d))
  753.     {
  754.       eneg (d);
  755.       s = 1;
  756.     }
  757.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  758.   ediv (df, d, dg);        /* dg = d / 2^32 is the high word */
  759.   euifrac (dg, high, dh);
  760.   emul (df, dh, dg);        /* fractional part is the low word */
  761.   euifrac (dg, low, dh);
  762.   if (s)
  763.     {
  764.       /* complement and add 1 */
  765.       *high = ~(*high);
  766.       if (*low)
  767.     *low = -(*low);
  768.       else
  769.     *high += 1;
  770.     }
  771. }
  772.  
  773.  
  774. /* REAL_VALUE_LDEXP macro.
  775.  */
  776. REAL_VALUE_TYPE
  777. ereal_ldexp (x, n)
  778.      REAL_VALUE_TYPE x;
  779.      int n;
  780. {
  781.   unsigned EMUSHORT e[NE], y[NE];
  782.   REAL_VALUE_TYPE r;
  783.  
  784.   GET_REAL (&x, e);
  785. #ifdef NANS
  786.   if (eisnan (e))
  787.     return (x);
  788. #endif
  789.   eldexp (e, n, y);
  790.   PUT_REAL (y, &r);
  791.   return (r);
  792. }
  793.  
  794. /* These routines are conditionally compiled because functions
  795.  * of the same names may be defined in fold-const.c.  */
  796. #ifdef REAL_ARITHMETIC
  797.  
  798. /* Check for infinity in a REAL_VALUE_TYPE. */
  799. int
  800. target_isinf (x)
  801.      REAL_VALUE_TYPE x;
  802. {
  803.   unsigned EMUSHORT e[NE];
  804.  
  805. #ifdef INFINITY
  806.   GET_REAL (&x, e);
  807.   return (eisinf (e));
  808. #else
  809.   return 0;
  810. #endif
  811. }
  812.  
  813.  
  814. /* Check whether a REAL_VALUE_TYPE item is a NaN. */
  815.  
  816. int
  817. target_isnan (x)
  818.      REAL_VALUE_TYPE x;
  819. {
  820.   unsigned EMUSHORT e[NE];
  821.  
  822. #ifdef NANS
  823.   GET_REAL (&x, e);
  824.   return (eisnan (e));
  825. #else
  826.   return (0);
  827. #endif
  828. }
  829.  
  830.  
  831. /* Check for a negative REAL_VALUE_TYPE number.
  832.  * this means strictly less than zero, not -0.
  833.  */
  834.  
  835. int
  836. target_negative (x)
  837.      REAL_VALUE_TYPE x;
  838. {
  839.   unsigned EMUSHORT e[NE];
  840.  
  841.   GET_REAL (&x, e);
  842.   if (ecmp (e, ezero) == -1)
  843.     return (1);
  844.   return (0);
  845. }
  846.  
  847. /* Expansion of REAL_VALUE_TRUNCATE.
  848.  * The result is in floating point, rounded to nearest or even.
  849.  */
  850. REAL_VALUE_TYPE
  851. real_value_truncate (mode, arg)
  852.      enum machine_mode mode;
  853.      REAL_VALUE_TYPE arg;
  854. {
  855.   unsigned EMUSHORT e[NE], t[NE];
  856.   REAL_VALUE_TYPE r;
  857.  
  858.   GET_REAL (&arg, e);
  859. #ifdef NANS
  860.   if (eisnan (e))
  861.     return (arg);
  862. #endif
  863.   eclear (t);
  864.   switch (mode)
  865.     {
  866.     case TFmode:
  867.       etoe113 (e, t);
  868.       e113toe (t, t);
  869.       break;
  870.  
  871.     case XFmode:
  872.       etoe64 (e, t);
  873.       e64toe (t, t);
  874.       break;
  875.  
  876.     case DFmode:
  877.       etoe53 (e, t);
  878.       e53toe (t, t);
  879.       break;
  880.  
  881.     case SFmode:
  882.       etoe24 (e, t);
  883.       e24toe (t, t);
  884.       break;
  885.  
  886.     case SImode:
  887.       r = etrunci (arg);
  888.       return (r);
  889.  
  890.     default:
  891.       abort ();
  892.     }
  893.   PUT_REAL (t, &r);
  894.   return (r);
  895. }
  896.  
  897. #endif /* REAL_ARITHMETIC defined */
  898.  
  899. /* Used for debugging--print the value of R in human-readable format
  900.    on stderr.  */
  901.  
  902. void
  903. debug_real (r)
  904.      REAL_VALUE_TYPE r;
  905. {
  906.   char dstr[30];
  907.  
  908.   REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
  909.   fprintf (stderr, "%s", dstr);
  910. }  
  911.  
  912.  
  913. /* Target values are arrays of host longs. A long is guaranteed
  914.    to be at least 32 bits wide. */
  915.  
  916. /* 128-bit long double */
  917. void 
  918. etartdouble (r, l)
  919.      REAL_VALUE_TYPE r;
  920.      long l[];
  921. {
  922.   unsigned EMUSHORT e[NE];
  923.  
  924.   GET_REAL (&r, e);
  925.   etoe113 (e, e);
  926.   endian (e, l, TFmode);
  927. }
  928.  
  929. /* 80-bit long double */
  930. void 
  931. etarldouble (r, l)
  932.      REAL_VALUE_TYPE r;
  933.      long l[];
  934. {
  935.   unsigned EMUSHORT e[NE];
  936.  
  937.   GET_REAL (&r, e);
  938.   etoe64 (e, e);
  939.   endian (e, l, XFmode);
  940. }
  941.  
  942. void 
  943. etardouble (r, l)
  944.      REAL_VALUE_TYPE r;
  945.      long l[];
  946. {
  947.   unsigned EMUSHORT e[NE];
  948.  
  949.   GET_REAL (&r, e);
  950.   etoe53 (e, e);
  951.   endian (e, l, DFmode);
  952. }
  953.  
  954. long
  955. etarsingle (r)
  956.      REAL_VALUE_TYPE r;
  957. {
  958.   unsigned EMUSHORT e[NE];
  959.   unsigned long l;
  960.  
  961.   GET_REAL (&r, e);
  962.   etoe24 (e, e);
  963.   endian (e, &l, SFmode);
  964.   return ((long) l);
  965. }
  966.  
  967. void
  968. ereal_to_decimal (x, s)
  969.      REAL_VALUE_TYPE x;
  970.      char *s;
  971. {
  972.   unsigned EMUSHORT e[NE];
  973.  
  974.   GET_REAL (&x, e);
  975.   etoasc (e, s, 20);
  976. }
  977.  
  978. int
  979. ereal_cmp (x, y)
  980.      REAL_VALUE_TYPE x, y;
  981. {
  982.   unsigned EMUSHORT ex[NE], ey[NE];
  983.  
  984.   GET_REAL (&x, ex);
  985.   GET_REAL (&y, ey);
  986.   return (ecmp (ex, ey));
  987. }
  988.  
  989. int
  990. ereal_isneg (x)
  991.      REAL_VALUE_TYPE x;
  992. {
  993.   unsigned EMUSHORT ex[NE];
  994.  
  995.   GET_REAL (&x, ex);
  996.   return (eisneg (ex));
  997. }
  998.  
  999. /* End of REAL_ARITHMETIC interface */
  1000.  
  1001. /*                            ieee.c
  1002.  *
  1003.  *    Extended precision IEEE binary floating point arithmetic routines
  1004.  *
  1005.  * Numbers are stored in C language as arrays of 16-bit unsigned
  1006.  * short integers.  The arguments of the routines are pointers to
  1007.  * the arrays.
  1008.  *
  1009.  *
  1010.  * External e type data structure, simulates Intel 8087 chip
  1011.  * temporary real format but possibly with a larger significand:
  1012.  *
  1013.  *    NE-1 significand words    (least significant word first,
  1014.  *                 most significant bit is normally set)
  1015.  *    exponent        (value = EXONE for 1.0,
  1016.  *                top bit is the sign)
  1017.  *
  1018.  *
  1019.  * Internal data structure of a number (a "word" is 16 bits):
  1020.  *
  1021.  * ei[0]    sign word    (0 for positive, 0xffff for negative)
  1022.  * ei[1]    biased exponent    (value = EXONE for the number 1.0)
  1023.  * ei[2]    high guard word    (always zero after normalization)
  1024.  * ei[3]
  1025.  * to ei[NI-2]    significand    (NI-4 significand words,
  1026.  *                 most significant word first,
  1027.  *                 most significant bit is set)
  1028.  * ei[NI-1]    low guard word    (0x8000 bit is rounding place)
  1029.  *
  1030.  *
  1031.  *
  1032.  *        Routines for external format numbers
  1033.  *
  1034.  *    asctoe (string, e)    ASCII string to extended double e type
  1035.  *    asctoe64 (string, &d)    ASCII string to long double
  1036.  *    asctoe53 (string, &d)    ASCII string to double
  1037.  *    asctoe24 (string, &f)    ASCII string to single
  1038.  *    asctoeg (string, e, prec) ASCII string to specified precision
  1039.  *    e24toe (&f, e)        IEEE single precision to e type
  1040.  *    e53toe (&d, e)        IEEE double precision to e type
  1041.  *    e64toe (&d, e)        IEEE long double precision to e type
  1042.  *    e113toe (&d, e)        128-bit long double precision to e type
  1043.  *    eabs (e)            absolute value
  1044.  *    eadd (a, b, c)        c = b + a
  1045.  *    eclear (e)        e = 0
  1046.  *    ecmp (a, b)        Returns 1 if a > b, 0 if a == b,
  1047.  *                -1 if a < b, -2 if either a or b is a NaN.
  1048.  *    ediv (a, b, c)        c = b / a
  1049.  *    efloor (a, b)        truncate to integer, toward -infinity
  1050.  *    efrexp (a, exp, s)    extract exponent and significand
  1051.  *    eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
  1052.  *    euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
  1053.  *    einfin (e)        set e to infinity, leaving its sign alone
  1054.  *    eldexp (a, n, b)    multiply by 2**n
  1055.  *    emov (a, b)        b = a
  1056.  *    emul (a, b, c)        c = b * a
  1057.  *    eneg (e)            e = -e
  1058.  *    eround (a, b)        b = nearest integer value to a
  1059.  *    esub (a, b, c)        c = b - a
  1060.  *    e24toasc (&f, str, n)    single to ASCII string, n digits after decimal
  1061.  *    e53toasc (&d, str, n)    double to ASCII string, n digits after decimal
  1062.  *    e64toasc (&d, str, n)    80-bit long double to ASCII string
  1063.  *    e113toasc (&d, str, n)    128-bit long double to ASCII string
  1064.  *    etoasc (e, str, n)    e to ASCII string, n digits after decimal
  1065.  *    etoe24 (e, &f)        convert e type to IEEE single precision
  1066.  *    etoe53 (e, &d)        convert e type to IEEE double precision
  1067.  *    etoe64 (e, &d)        convert e type to IEEE long double precision
  1068.  *    ltoe (&l, e)        HOST_WIDE_INT to e type
  1069.  *    ultoe (&l, e)        unsigned HOST_WIDE_INT to e type
  1070.  *      eisneg (e)              1 if sign bit of e != 0, else 0
  1071.  *      eisinf (e)              1 if e has maximum exponent (non-IEEE)
  1072.  *                or is infinite (IEEE)
  1073.  *      eisnan (e)              1 if e is a NaN
  1074.  *
  1075.  *
  1076.  *        Routines for internal format numbers
  1077.  *
  1078.  *    eaddm (ai, bi)        add significands, bi = bi + ai
  1079.  *    ecleaz (ei)        ei = 0
  1080.  *    ecleazs (ei)        set ei = 0 but leave its sign alone
  1081.  *    ecmpm (ai, bi)        compare significands, return 1, 0, or -1
  1082.  *    edivm (ai, bi)        divide  significands, bi = bi / ai
  1083.  *    emdnorm (ai,l,s,exp)    normalize and round off
  1084.  *    emovi (a, ai)        convert external a to internal ai
  1085.  *    emovo (ai, a)        convert internal ai to external a
  1086.  *    emovz (ai, bi)        bi = ai, low guard word of bi = 0
  1087.  *    emulm (ai, bi)        multiply significands, bi = bi * ai
  1088.  *    enormlz (ei)        left-justify the significand
  1089.  *    eshdn1 (ai)        shift significand and guards down 1 bit
  1090.  *    eshdn8 (ai)        shift down 8 bits
  1091.  *    eshdn6 (ai)        shift down 16 bits
  1092.  *    eshift (ai, n)        shift ai n bits up (or down if n < 0)
  1093.  *    eshup1 (ai)        shift significand and guards up 1 bit
  1094.  *    eshup8 (ai)        shift up 8 bits
  1095.  *    eshup6 (ai)        shift up 16 bits
  1096.  *    esubm (ai, bi)        subtract significands, bi = bi - ai
  1097.  *      eiisinf (ai)            1 if infinite
  1098.  *      eiisnan (ai)            1 if a NaN
  1099.  *      einan (ai)              set ai = NaN
  1100.  *      eiinfin (ai)            set ai = infinity
  1101.  *
  1102.  *
  1103.  * The result is always normalized and rounded to NI-4 word precision
  1104.  * after each arithmetic operation.
  1105.  *
  1106.  * Exception flags are NOT fully supported.
  1107.  *
  1108.  * Signaling NaN's are NOT supported; they are treated the same
  1109.  * as quiet NaN's.
  1110.  *
  1111.  * Define INFINITY for support of infinity; otherwise a
  1112.  * saturation arithmetic is implemented.
  1113.  *
  1114.  * Define NANS for support of Not-a-Number items; otherwise the
  1115.  * arithmetic will never produce a NaN output, and might be confused
  1116.  * by a NaN input.
  1117.  * If NaN's are supported, the output of `ecmp (a,b)' is -2 if
  1118.  * either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
  1119.  * may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
  1120.  * if in doubt.
  1121.  *
  1122.  * Denormals are always supported here where appropriate (e.g., not
  1123.  * for conversion to DEC numbers).
  1124.  *
  1125.  */
  1126.  
  1127.  
  1128. /*                            mconf.h
  1129.  *
  1130.  *    Common include file for math routines
  1131.  *
  1132.  *
  1133.  *
  1134.  * SYNOPSIS:
  1135.  *
  1136.  * #include "mconf.h"
  1137.  *
  1138.  *
  1139.  *
  1140.  * DESCRIPTION:
  1141.  *
  1142.  * This file contains definitions for error codes that are
  1143.  * passed to the common error handling routine mtherr
  1144.  * (which see).
  1145.  *
  1146.  * The file also includes a conditional assembly definition
  1147.  * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
  1148.  * IEEE, or UNKnown).
  1149.  *
  1150.  * For Digital Equipment PDP-11 and VAX computers, certain
  1151.  * IBM systems, and others that use numbers with a 56-bit
  1152.  * significand, the symbol DEC should be defined.  In this
  1153.  * mode, most floating point constants are given as arrays
  1154.  * of octal integers to eliminate decimal to binary conversion
  1155.  * errors that might be introduced by the compiler.
  1156.  *
  1157.  * For computers, such as IBM PC, that follow the IEEE
  1158.  * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
  1159.  * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
  1160.  * These numbers have 53-bit significands.  In this mode, constants
  1161.  * are provided as arrays of hexadecimal 16 bit integers.
  1162.  *
  1163.  * To accommodate other types of computer arithmetic, all
  1164.  * constants are also provided in a normal decimal radix
  1165.  * which one can hope are correctly converted to a suitable
  1166.  * format by the available C language compiler.  To invoke
  1167.  * this mode, the symbol UNK is defined.
  1168.  *
  1169.  * An important difference among these modes is a predefined
  1170.  * set of machine arithmetic constants for each.  The numbers
  1171.  * MACHEP (the machine roundoff error), MAXNUM (largest number
  1172.  * represented), and several other parameters are preset by
  1173.  * the configuration symbol.  Check the file const.c to
  1174.  * ensure that these values are correct for your computer.
  1175.  *
  1176.  * For ANSI C compatibility, define ANSIC equal to 1.  Currently
  1177.  * this affects only the atan2 function and others that use it.
  1178.  */
  1179.  
  1180. /* Constant definitions for math error conditions.  */
  1181.  
  1182. #define DOMAIN        1    /* argument domain error */
  1183. #define SING        2    /* argument singularity */
  1184. #define OVERFLOW    3    /* overflow range error */
  1185. #define UNDERFLOW    4    /* underflow range error */
  1186. #define TLOSS        5    /* total loss of precision */
  1187. #define PLOSS        6    /* partial loss of precision */
  1188. #define INVALID        7    /* NaN-producing operation */
  1189.  
  1190. /*  e type constants used by high precision check routines */
  1191.  
  1192. #if LONG_DOUBLE_TYPE_SIZE == 128
  1193. /* 0.0 */
  1194. unsigned EMUSHORT ezero[NE] =
  1195.  {0x0000, 0x0000, 0x0000, 0x0000,
  1196.   0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
  1197. extern unsigned EMUSHORT ezero[];
  1198.  
  1199. /* 5.0E-1 */
  1200. unsigned EMUSHORT ehalf[NE] =
  1201.  {0x0000, 0x0000, 0x0000, 0x0000,
  1202.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
  1203. extern unsigned EMUSHORT ehalf[];
  1204.  
  1205. /* 1.0E0 */
  1206. unsigned EMUSHORT eone[NE] =
  1207.  {0x0000, 0x0000, 0x0000, 0x0000,
  1208.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
  1209. extern unsigned EMUSHORT eone[];
  1210.  
  1211. /* 2.0E0 */
  1212. unsigned EMUSHORT etwo[NE] =
  1213.  {0x0000, 0x0000, 0x0000, 0x0000,
  1214.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
  1215. extern unsigned EMUSHORT etwo[];
  1216.  
  1217. /* 3.2E1 */
  1218. unsigned EMUSHORT e32[NE] =
  1219.  {0x0000, 0x0000, 0x0000, 0x0000,
  1220.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
  1221. extern unsigned EMUSHORT e32[];
  1222.  
  1223. /* 6.93147180559945309417232121458176568075500134360255E-1 */
  1224. unsigned EMUSHORT elog2[NE] =
  1225.  {0x40f3, 0xf6af, 0x03f2, 0xb398,
  1226.   0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
  1227. extern unsigned EMUSHORT elog2[];
  1228.  
  1229. /* 1.41421356237309504880168872420969807856967187537695E0 */
  1230. unsigned EMUSHORT esqrt2[NE] =
  1231.  {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
  1232.   0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
  1233. extern unsigned EMUSHORT esqrt2[];
  1234.  
  1235. /* 3.14159265358979323846264338327950288419716939937511E0 */
  1236. unsigned EMUSHORT epi[NE] =
  1237.  {0x2902, 0x1cd1, 0x80dc, 0x628b,
  1238.   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
  1239. extern unsigned EMUSHORT epi[];
  1240.  
  1241. #else
  1242. /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
  1243. unsigned EMUSHORT ezero[NE] =
  1244.  {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
  1245. unsigned EMUSHORT ehalf[NE] =
  1246.  {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
  1247. unsigned EMUSHORT eone[NE] =
  1248.  {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
  1249. unsigned EMUSHORT etwo[NE] =
  1250.  {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
  1251. unsigned EMUSHORT e32[NE] =
  1252.  {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
  1253. unsigned EMUSHORT elog2[NE] =
  1254.  {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
  1255. unsigned EMUSHORT esqrt2[NE] =
  1256.  {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
  1257. unsigned EMUSHORT epi[NE] =
  1258.  {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
  1259. #endif
  1260.  
  1261.  
  1262.  
  1263. /* Control register for rounding precision.
  1264.  * This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.
  1265.  */
  1266. int rndprc = NBITS;
  1267. extern int rndprc;
  1268.  
  1269. void eaddm (), esubm (), emdnorm (), asctoeg ();
  1270. static void toe24 (), toe53 (), toe64 (), toe113 ();
  1271. void eremain (), einit (), eiremain ();
  1272. int ecmpm (), edivm (), emulm ();
  1273. void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
  1274. #ifdef DEC
  1275. void etodec (), todec (), dectoe ();
  1276. #endif
  1277. #ifdef IBM
  1278. void etoibm (), toibm (), ibmtoe ();
  1279. #endif
  1280.  
  1281.  
  1282. void 
  1283. einit ()
  1284. {
  1285. }
  1286.  
  1287. /*
  1288. ; Clear out entire external format number.
  1289. ;
  1290. ; unsigned EMUSHORT x[];
  1291. ; eclear (x);
  1292. */
  1293.  
  1294. void 
  1295. eclear (x)
  1296.      register unsigned EMUSHORT *x;
  1297. {
  1298.   register int i;
  1299.  
  1300.   for (i = 0; i < NE; i++)
  1301.     *x++ = 0;
  1302. }
  1303.  
  1304.  
  1305.  
  1306. /* Move external format number from a to b.
  1307.  *
  1308.  * emov (a, b);
  1309.  */
  1310.  
  1311. void 
  1312. emov (a, b)
  1313.      register unsigned EMUSHORT *a, *b;
  1314. {
  1315.   register int i;
  1316.  
  1317.   for (i = 0; i < NE; i++)
  1318.     *b++ = *a++;
  1319. }
  1320.  
  1321.  
  1322. /*
  1323. ;    Absolute value of external format number
  1324. ;
  1325. ;    EMUSHORT x[NE];
  1326. ;    eabs (x);
  1327. */
  1328.  
  1329. void 
  1330. eabs (x)
  1331.      unsigned EMUSHORT x[];    /* x is the memory address of a short */
  1332. {
  1333.  
  1334.   x[NE - 1] &= 0x7fff;        /* sign is top bit of last word of external format */
  1335. }
  1336.  
  1337.  
  1338.  
  1339.  
  1340. /*
  1341. ;    Negate external format number
  1342. ;
  1343. ;    unsigned EMUSHORT x[NE];
  1344. ;    eneg (x);
  1345. */
  1346.  
  1347. void 
  1348. eneg (x)
  1349.      unsigned EMUSHORT x[];
  1350. {
  1351.  
  1352. #ifdef NANS
  1353.   if (eisnan (x))
  1354.     return;
  1355. #endif
  1356.   x[NE - 1] ^= 0x8000;        /* Toggle the sign bit */
  1357. }
  1358.  
  1359.  
  1360.  
  1361. /* Return 1 if external format number is negative,
  1362.  * else return zero, including when it is a NaN.
  1363.  */
  1364. int 
  1365. eisneg (x)
  1366.      unsigned EMUSHORT x[];
  1367. {
  1368.  
  1369. #ifdef NANS
  1370.   if (eisnan (x))
  1371.     return (0);
  1372. #endif
  1373.   if (x[NE - 1] & 0x8000)
  1374.     return (1);
  1375.   else
  1376.     return (0);
  1377. }
  1378.  
  1379.  
  1380. /* Return 1 if external format number is infinity.
  1381.  * else return zero.
  1382.  */
  1383. int 
  1384. eisinf (x)
  1385.      unsigned EMUSHORT x[];
  1386. {
  1387.  
  1388. #ifdef NANS
  1389.   if (eisnan (x))
  1390.     return (0);
  1391. #endif
  1392.   if ((x[NE - 1] & 0x7fff) == 0x7fff)
  1393.     return (1);
  1394.   else
  1395.     return (0);
  1396. }
  1397.  
  1398.  
  1399. /* Check if e-type number is not a number.
  1400.    The bit pattern is one that we defined, so we know for sure how to
  1401.    detect it.  */
  1402.  
  1403. int 
  1404. eisnan (x)
  1405.      unsigned EMUSHORT x[];
  1406. {
  1407.  
  1408. #ifdef NANS
  1409.   int i;
  1410. /* NaN has maximum exponent */
  1411.   if ((x[NE - 1] & 0x7fff) != 0x7fff)
  1412.     return (0);
  1413. /* ... and non-zero significand field. */
  1414.   for (i = 0; i < NE - 1; i++)
  1415.     {
  1416.       if (*x++ != 0)
  1417.         return (1);
  1418.     }
  1419. #endif
  1420.   return (0);
  1421. }
  1422.  
  1423. /*  Fill external format number with infinity pattern (IEEE)
  1424.     or largest possible number (non-IEEE). */
  1425.  
  1426. void 
  1427. einfin (x)
  1428.      register unsigned EMUSHORT *x;
  1429. {
  1430.   register int i;
  1431.  
  1432. #ifdef INFINITY
  1433.   for (i = 0; i < NE - 1; i++)
  1434.     *x++ = 0;
  1435.   *x |= 32767;
  1436. #else
  1437.   for (i = 0; i < NE - 1; i++)
  1438.     *x++ = 0xffff;
  1439.   *x |= 32766;
  1440.   if (rndprc < NBITS)
  1441.     {
  1442.       if (rndprc == 113)
  1443.     {
  1444.       *(x - 9) = 0;
  1445.       *(x - 8) = 0;
  1446.     }
  1447.       if (rndprc == 64)
  1448.     {
  1449.       *(x - 5) = 0;
  1450.     }
  1451.       if (rndprc == 53)
  1452.     {
  1453.       *(x - 4) = 0xf800;
  1454.     }
  1455.       else
  1456.     {
  1457.       *(x - 4) = 0;
  1458.       *(x - 3) = 0;
  1459.       *(x - 2) = 0xff00;
  1460.     }
  1461.     }
  1462. #endif
  1463. }
  1464.  
  1465.  
  1466. /* Output an e-type NaN.
  1467.    This generates Intel's quiet NaN pattern for extended real.
  1468.    The exponent is 7fff, the leading mantissa word is c000.  */
  1469.  
  1470. void 
  1471. enan (x)
  1472.      register unsigned EMUSHORT *x;
  1473. {
  1474.   register int i;
  1475.  
  1476.   for (i = 0; i < NE - 2; i++)
  1477.     *x++ = 0;
  1478.   *x++ = 0xc000;
  1479.   *x = 0x7fff;
  1480. }
  1481.  
  1482.  
  1483. /* Move in external format number,
  1484.  * converting it to internal format.
  1485.  */
  1486. void 
  1487. emovi (a, b)
  1488.      unsigned EMUSHORT *a, *b;
  1489. {
  1490.   register unsigned EMUSHORT *p, *q;
  1491.   int i;
  1492.  
  1493.   q = b;
  1494.   p = a + (NE - 1);        /* point to last word of external number */
  1495.   /* get the sign bit */
  1496.   if (*p & 0x8000)
  1497.     *q++ = 0xffff;
  1498.   else
  1499.     *q++ = 0;
  1500.   /* get the exponent */
  1501.   *q = *p--;
  1502.   *q++ &= 0x7fff;        /* delete the sign bit */
  1503. #ifdef INFINITY
  1504.   if ((*(q - 1) & 0x7fff) == 0x7fff)
  1505.     {
  1506. #ifdef NANS
  1507.       if (eisnan (a))
  1508.     {
  1509.       *q++ = 0;
  1510.       for (i = 3; i < NI; i++)
  1511.         *q++ = *p--;
  1512.       return;
  1513.     }
  1514. #endif
  1515.       for (i = 2; i < NI; i++)
  1516.     *q++ = 0;
  1517.       return;
  1518.     }
  1519. #endif
  1520.   /* clear high guard word */
  1521.   *q++ = 0;
  1522.   /* move in the significand */
  1523.   for (i = 0; i < NE - 1; i++)
  1524.     *q++ = *p--;
  1525.   /* clear low guard word */
  1526.   *q = 0;
  1527. }
  1528.  
  1529.  
  1530. /* Move internal format number out,
  1531.  * converting it to external format.
  1532.  */
  1533. void 
  1534. emovo (a, b)
  1535.      unsigned EMUSHORT *a, *b;
  1536. {
  1537.   register unsigned EMUSHORT *p, *q;
  1538.   unsigned EMUSHORT i;
  1539.  
  1540.   p = a;
  1541.   q = b + (NE - 1);        /* point to output exponent */
  1542.   /* combine sign and exponent */
  1543.   i = *p++;
  1544.   if (i)
  1545.     *q-- = *p++ | 0x8000;
  1546.   else
  1547.     *q-- = *p++;
  1548. #ifdef INFINITY
  1549.   if (*(p - 1) == 0x7fff)
  1550.     {
  1551. #ifdef NANS
  1552.       if (eiisnan (a))
  1553.     {
  1554.       enan (b);
  1555.       return;
  1556.     }
  1557. #endif
  1558.       einfin (b);
  1559.     return;
  1560.     }
  1561. #endif
  1562.   /* skip over guard word */
  1563.   ++p;
  1564.   /* move the significand */
  1565.   for (i = 0; i < NE - 1; i++)
  1566.     *q-- = *p++;
  1567. }
  1568.  
  1569.  
  1570.  
  1571.  
  1572. /* Clear out internal format number.
  1573.  */
  1574.  
  1575. void 
  1576. ecleaz (xi)
  1577.      register unsigned EMUSHORT *xi;
  1578. {
  1579.   register int i;
  1580.  
  1581.   for (i = 0; i < NI; i++)
  1582.     *xi++ = 0;
  1583. }
  1584.  
  1585.  
  1586. /* same, but don't touch the sign. */
  1587.  
  1588. void 
  1589. ecleazs (xi)
  1590.      register unsigned EMUSHORT *xi;
  1591. {
  1592.   register int i;
  1593.  
  1594.   ++xi;
  1595.   for (i = 0; i < NI - 1; i++)
  1596.     *xi++ = 0;
  1597. }
  1598.  
  1599.  
  1600.  
  1601. /* Move internal format number from a to b.
  1602.  */
  1603. void 
  1604. emovz (a, b)
  1605.      register unsigned EMUSHORT *a, *b;
  1606. {
  1607.   register int i;
  1608.  
  1609.   for (i = 0; i < NI - 1; i++)
  1610.     *b++ = *a++;
  1611.   /* clear low guard word */
  1612.   *b = 0;
  1613. }
  1614.  
  1615. /* Generate internal format NaN.
  1616.    The explicit pattern for this is maximum exponent and
  1617.    top two significand bits set.  */
  1618.  
  1619. void
  1620. einan (x)
  1621.      unsigned EMUSHORT x[];
  1622. {
  1623.  
  1624.   ecleaz (x);
  1625.   x[E] = 0x7fff;
  1626.   x[M + 1] = 0xc000;
  1627. }
  1628.  
  1629. /* Return nonzero if internal format number is a NaN. */
  1630.  
  1631. int 
  1632. eiisnan (x)
  1633.      unsigned EMUSHORT x[];
  1634. {
  1635.   int i;
  1636.  
  1637.   if ((x[E] & 0x7fff) == 0x7fff)
  1638.     {
  1639.       for (i = M + 1; i < NI; i++)
  1640.     {
  1641.       if (x[i] != 0)
  1642.         return (1);
  1643.     }
  1644.     }
  1645.   return (0);
  1646. }
  1647.  
  1648. /* Fill internal format number with infinity pattern.
  1649.    This has maximum exponent and significand all zeros.  */
  1650.  
  1651. void
  1652. eiinfin (x)
  1653.      unsigned EMUSHORT x[];
  1654. {
  1655.  
  1656.   ecleaz (x);
  1657.   x[E] = 0x7fff;
  1658. }
  1659.  
  1660. /* Return nonzero if internal format number is infinite. */
  1661.  
  1662. int 
  1663. eiisinf (x)
  1664.      unsigned EMUSHORT x[];
  1665. {
  1666.  
  1667. #ifdef NANS
  1668.   if (eiisnan (x))
  1669.     return (0);
  1670. #endif
  1671.   if ((x[E] & 0x7fff) == 0x7fff)
  1672.     return (1);
  1673.   return (0);
  1674. }
  1675.  
  1676.  
  1677. /*
  1678. ;    Compare significands of numbers in internal format.
  1679. ;    Guard words are included in the comparison.
  1680. ;
  1681. ;    unsigned EMUSHORT a[NI], b[NI];
  1682. ;    cmpm (a, b);
  1683. ;
  1684. ;    for the significands:
  1685. ;    returns    +1 if a > b
  1686. ;         0 if a == b
  1687. ;        -1 if a < b
  1688. */
  1689. int
  1690. ecmpm (a, b)
  1691.      register unsigned EMUSHORT *a, *b;
  1692. {
  1693.   int i;
  1694.  
  1695.   a += M;            /* skip up to significand area */
  1696.   b += M;
  1697.   for (i = M; i < NI; i++)
  1698.     {
  1699.       if (*a++ != *b++)
  1700.     goto difrnt;
  1701.     }
  1702.   return (0);
  1703.  
  1704.  difrnt:
  1705.   if (*(--a) > *(--b))
  1706.     return (1);
  1707.   else
  1708.     return (-1);
  1709. }
  1710.  
  1711.  
  1712. /*
  1713. ;    Shift significand down by 1 bit
  1714. */
  1715.  
  1716. void 
  1717. eshdn1 (x)
  1718.      register unsigned EMUSHORT *x;
  1719. {
  1720.   register unsigned EMUSHORT bits;
  1721.   int i;
  1722.  
  1723.   x += M;            /* point to significand area */
  1724.  
  1725.   bits = 0;
  1726.   for (i = M; i < NI; i++)
  1727.     {
  1728.       if (*x & 1)
  1729.     bits |= 1;
  1730.       *x >>= 1;
  1731.       if (bits & 2)
  1732.     *x |= 0x8000;
  1733.       bits <<= 1;
  1734.       ++x;
  1735.     }
  1736. }
  1737.  
  1738.  
  1739.  
  1740. /*
  1741. ;    Shift significand up by 1 bit
  1742. */
  1743.  
  1744. void 
  1745. eshup1 (x)
  1746.      register unsigned EMUSHORT *x;
  1747. {
  1748.   register unsigned EMUSHORT bits;
  1749.   int i;
  1750.  
  1751.   x += NI - 1;
  1752.   bits = 0;
  1753.  
  1754.   for (i = M; i < NI; i++)
  1755.     {
  1756.       if (*x & 0x8000)
  1757.     bits |= 1;
  1758.       *x <<= 1;
  1759.       if (bits & 2)
  1760.     *x |= 1;
  1761.       bits <<= 1;
  1762.       --x;
  1763.     }
  1764. }
  1765.  
  1766.  
  1767.  
  1768. /*
  1769. ;    Shift significand down by 8 bits
  1770. */
  1771.  
  1772. void 
  1773. eshdn8 (x)
  1774.      register unsigned EMUSHORT *x;
  1775. {
  1776.   register unsigned EMUSHORT newbyt, oldbyt;
  1777.   int i;
  1778.  
  1779.   x += M;
  1780.   oldbyt = 0;
  1781.   for (i = M; i < NI; i++)
  1782.     {
  1783.       newbyt = *x << 8;
  1784.       *x >>= 8;
  1785.       *x |= oldbyt;
  1786.       oldbyt = newbyt;
  1787.       ++x;
  1788.     }
  1789. }
  1790.  
  1791. /*
  1792. ;    Shift significand up by 8 bits
  1793. */
  1794.  
  1795. void 
  1796. eshup8 (x)
  1797.      register unsigned EMUSHORT *x;
  1798. {
  1799.   int i;
  1800.   register unsigned EMUSHORT newbyt, oldbyt;
  1801.  
  1802.   x += NI - 1;
  1803.   oldbyt = 0;
  1804.  
  1805.   for (i = M; i < NI; i++)
  1806.     {
  1807.       newbyt = *x >> 8;
  1808.       *x <<= 8;
  1809.       *x |= oldbyt;
  1810.       oldbyt = newbyt;
  1811.       --x;
  1812.     }
  1813. }
  1814.  
  1815. /*
  1816. ;    Shift significand up by 16 bits
  1817. */
  1818.  
  1819. void 
  1820. eshup6 (x)
  1821.      register unsigned EMUSHORT *x;
  1822. {
  1823.   int i;
  1824.   register unsigned EMUSHORT *p;
  1825.  
  1826.   p = x + M;
  1827.   x += M + 1;
  1828.  
  1829.   for (i = M; i < NI - 1; i++)
  1830.     *p++ = *x++;
  1831.  
  1832.   *p = 0;
  1833. }
  1834.  
  1835. /*
  1836. ;    Shift significand down by 16 bits
  1837. */
  1838.  
  1839. void 
  1840. eshdn6 (x)
  1841.      register unsigned EMUSHORT *x;
  1842. {
  1843.   int i;
  1844.   register unsigned EMUSHORT *p;
  1845.  
  1846.   x += NI - 1;
  1847.   p = x + 1;
  1848.  
  1849.   for (i = M; i < NI - 1; i++)
  1850.     *(--p) = *(--x);
  1851.  
  1852.   *(--p) = 0;
  1853. }
  1854.  
  1855. /*
  1856. ;    Add significands
  1857. ;    x + y replaces y
  1858. */
  1859.  
  1860. void 
  1861. eaddm (x, y)
  1862.      unsigned EMUSHORT *x, *y;
  1863. {
  1864.   register unsigned EMULONG a;
  1865.   int i;
  1866.   unsigned int carry;
  1867.  
  1868.   x += NI - 1;
  1869.   y += NI - 1;
  1870.   carry = 0;
  1871.   for (i = M; i < NI; i++)
  1872.     {
  1873.       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
  1874.       if (a & 0x10000)
  1875.     carry = 1;
  1876.       else
  1877.     carry = 0;
  1878.       *y = (unsigned EMUSHORT) a;
  1879.       --x;
  1880.       --y;
  1881.     }
  1882. }
  1883.  
  1884. /*
  1885. ;    Subtract significands
  1886. ;    y - x replaces y
  1887. */
  1888.  
  1889. void 
  1890. esubm (x, y)
  1891.      unsigned EMUSHORT *x, *y;
  1892. {
  1893.   unsigned EMULONG a;
  1894.   int i;
  1895.   unsigned int carry;
  1896.  
  1897.   x += NI - 1;
  1898.   y += NI - 1;
  1899.   carry = 0;
  1900.   for (i = M; i < NI; i++)
  1901.     {
  1902.       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
  1903.       if (a & 0x10000)
  1904.     carry = 1;
  1905.       else
  1906.     carry = 0;
  1907.       *y = (unsigned EMUSHORT) a;
  1908.       --x;
  1909.       --y;
  1910.     }
  1911. }
  1912.  
  1913.  
  1914. static unsigned EMUSHORT equot[NI];
  1915.  
  1916.  
  1917. #if 0
  1918. /* Radix 2 shift-and-add versions of multiply and divide  */
  1919.  
  1920.  
  1921. /* Divide significands */
  1922.  
  1923. int 
  1924. edivm (den, num)
  1925.      unsigned EMUSHORT den[], num[];
  1926. {
  1927.   int i;
  1928.   register unsigned EMUSHORT *p, *q;
  1929.   unsigned EMUSHORT j;
  1930.  
  1931.   p = &equot[0];
  1932.   *p++ = num[0];
  1933.   *p++ = num[1];
  1934.  
  1935.   for (i = M; i < NI; i++)
  1936.     {
  1937.       *p++ = 0;
  1938.     }
  1939.  
  1940.   /* Use faster compare and subtraction if denominator
  1941.    * has only 15 bits of significance.
  1942.    */
  1943.   p = &den[M + 2];
  1944.   if (*p++ == 0)
  1945.     {
  1946.       for (i = M + 3; i < NI; i++)
  1947.     {
  1948.       if (*p++ != 0)
  1949.         goto fulldiv;
  1950.     }
  1951.       if ((den[M + 1] & 1) != 0)
  1952.     goto fulldiv;
  1953.       eshdn1 (num);
  1954.       eshdn1 (den);
  1955.  
  1956.       p = &den[M + 1];
  1957.       q = &num[M + 1];
  1958.  
  1959.       for (i = 0; i < NBITS + 2; i++)
  1960.     {
  1961.       if (*p <= *q)
  1962.         {
  1963.           *q -= *p;
  1964.           j = 1;
  1965.         }
  1966.       else
  1967.         {
  1968.           j = 0;
  1969.         }
  1970.       eshup1 (equot);
  1971.       equot[NI - 2] |= j;
  1972.       eshup1 (num);
  1973.     }
  1974.       goto divdon;
  1975.     }
  1976.  
  1977.   /* The number of quotient bits to calculate is
  1978.    * NBITS + 1 scaling guard bit + 1 roundoff bit.
  1979.    */
  1980.  fulldiv:
  1981.  
  1982.   p = &equot[NI - 2];
  1983.   for (i = 0; i < NBITS + 2; i++)
  1984.     {
  1985.       if (ecmpm (den, num) <= 0)
  1986.     {
  1987.       esubm (den, num);
  1988.       j = 1;        /* quotient bit = 1 */
  1989.     }
  1990.       else
  1991.     j = 0;
  1992.       eshup1 (equot);
  1993.       *p |= j;
  1994.       eshup1 (num);
  1995.     }
  1996.  
  1997.  divdon:
  1998.  
  1999.   eshdn1 (equot);
  2000.   eshdn1 (equot);
  2001.  
  2002.   /* test for nonzero remainder after roundoff bit */
  2003.   p = &num[M];
  2004.   j = 0;
  2005.   for (i = M; i < NI; i++)
  2006.     {
  2007.       j |= *p++;
  2008.     }
  2009.   if (j)
  2010.     j = 1;
  2011.  
  2012.  
  2013.   for (i = 0; i < NI; i++)
  2014.     num[i] = equot[i];
  2015.   return ((int) j);
  2016. }
  2017.  
  2018.  
  2019. /* Multiply significands */
  2020. int 
  2021. emulm (a, b)
  2022.      unsigned EMUSHORT a[], b[];
  2023. {
  2024.   unsigned EMUSHORT *p, *q;
  2025.   int i, j, k;
  2026.  
  2027.   equot[0] = b[0];
  2028.   equot[1] = b[1];
  2029.   for (i = M; i < NI; i++)
  2030.     equot[i] = 0;
  2031.  
  2032.   p = &a[NI - 2];
  2033.   k = NBITS;
  2034.   while (*p == 0)        /* significand is not supposed to be all zero */
  2035.     {
  2036.       eshdn6 (a);
  2037.       k -= 16;
  2038.     }
  2039.   if ((*p & 0xff) == 0)
  2040.     {
  2041.       eshdn8 (a);
  2042.       k -= 8;
  2043.     }
  2044.  
  2045.   q = &equot[NI - 1];
  2046.   j = 0;
  2047.   for (i = 0; i < k; i++)
  2048.     {
  2049.       if (*p & 1)
  2050.     eaddm (b, equot);
  2051.       /* remember if there were any nonzero bits shifted out */
  2052.       if (*q & 1)
  2053.     j |= 1;
  2054.       eshdn1 (a);
  2055.       eshdn1 (equot);
  2056.     }
  2057.  
  2058.   for (i = 0; i < NI; i++)
  2059.     b[i] = equot[i];
  2060.  
  2061.   /* return flag for lost nonzero bits */
  2062.   return (j);
  2063. }
  2064.  
  2065. #else
  2066.  
  2067. /* Radix 65536 versions of multiply and divide  */
  2068.  
  2069.  
  2070. /* Multiply significand of e-type number b
  2071. by 16-bit quantity a, e-type result to c. */
  2072.  
  2073. void
  2074. m16m (a, b, c)
  2075.      unsigned short a;
  2076.      unsigned short b[], c[];
  2077. {
  2078.   register unsigned short *pp;
  2079.   register unsigned long carry;
  2080.   unsigned short *ps;
  2081.   unsigned short p[NI];
  2082.   unsigned long aa, m;
  2083.   int i;
  2084.  
  2085.   aa = a;
  2086.   pp = &p[NI-2];
  2087.   *pp++ = 0;
  2088.   *pp = 0;
  2089.   ps = &b[NI-1];
  2090.  
  2091.   for (i=M+1; i<NI; i++)
  2092.     {
  2093.       if (*ps == 0)
  2094.     {
  2095.       --ps;
  2096.       --pp;
  2097.       *(pp-1) = 0;
  2098.     }
  2099.       else
  2100.     {
  2101.       m = (unsigned long) aa * *ps--;
  2102.       carry = (m & 0xffff) + *pp;
  2103.       *pp-- = (unsigned short)carry;
  2104.       carry = (carry >> 16) + (m >> 16) + *pp;
  2105.       *pp = (unsigned short)carry;
  2106.       *(pp-1) = carry >> 16;
  2107.     }
  2108.     }
  2109.   for (i=M; i<NI; i++)
  2110.     c[i] = p[i];
  2111. }
  2112.  
  2113.  
  2114. /* Divide significands. Neither the numerator nor the denominator
  2115.    is permitted to have its high guard word nonzero.  */
  2116.  
  2117. int
  2118. edivm (den, num)
  2119.      unsigned short den[], num[];
  2120. {
  2121.   int i;
  2122.   register unsigned short *p;
  2123.   unsigned long tnum;
  2124.   unsigned short j, tdenm, tquot;
  2125.   unsigned short tprod[NI+1];
  2126.  
  2127.   p = &equot[0];
  2128.   *p++ = num[0];
  2129.   *p++ = num[1];
  2130.  
  2131.   for (i=M; i<NI; i++)
  2132.     {
  2133.       *p++ = 0;
  2134.     }
  2135.   eshdn1 (num);
  2136.   tdenm = den[M+1];
  2137.   for (i=M; i<NI; i++)
  2138.     {
  2139.       /* Find trial quotient digit (the radix is 65536). */
  2140.       tnum = (((unsigned long) num[M]) << 16) + num[M+1];
  2141.  
  2142.       /* Do not execute the divide instruction if it will overflow. */
  2143.       if ((tdenm * 0xffffL) < tnum)
  2144.     tquot = 0xffff;
  2145.       else
  2146.     tquot = tnum / tdenm;
  2147.       /* Multiply denominator by trial quotient digit. */
  2148.       m16m (tquot, den, tprod);
  2149.       /* The quotient digit may have been overestimated. */
  2150.       if (ecmpm (tprod, num) > 0)
  2151.     {
  2152.       tquot -= 1;
  2153.       esubm (den, tprod);
  2154.       if (ecmpm (tprod, num) > 0)
  2155.         {
  2156.           tquot -= 1;
  2157.           esubm (den, tprod);
  2158.         }
  2159.     }
  2160.       esubm (tprod, num);
  2161.       equot[i] = tquot;
  2162.       eshup6(num);
  2163.     }
  2164.   /* test for nonzero remainder after roundoff bit */
  2165.   p = &num[M];
  2166.   j = 0;
  2167.   for (i=M; i<NI; i++)
  2168.     {
  2169.       j |= *p++;
  2170.     }
  2171.   if (j)
  2172.     j = 1;
  2173.  
  2174.   for (i=0; i<NI; i++)
  2175.     num[i] = equot[i];
  2176.  
  2177.   return ((int)j);
  2178. }
  2179.  
  2180.  
  2181.  
  2182. /* Multiply significands */
  2183. int
  2184. emulm (a, b)
  2185.      unsigned short a[], b[];
  2186. {
  2187.   unsigned short *p, *q;
  2188.   unsigned short pprod[NI];
  2189.   unsigned short j;
  2190.   int i;
  2191.  
  2192.   equot[0] = b[0];
  2193.   equot[1] = b[1];
  2194.   for (i=M; i<NI; i++)
  2195.     equot[i] = 0;
  2196.  
  2197.   j = 0;
  2198.   p = &a[NI-1];
  2199.   q = &equot[NI-1];
  2200.   for (i=M+1; i<NI; i++)
  2201.     {
  2202.       if (*p == 0)
  2203.     {
  2204.       --p;
  2205.     }
  2206.       else
  2207.     {
  2208.       m16m (*p--, b, pprod);
  2209.       eaddm(pprod, equot);
  2210.     }
  2211.       j |= *q;
  2212.       eshdn6(equot);
  2213.     }
  2214.  
  2215.   for (i=0; i<NI; i++)
  2216.     b[i] = equot[i];
  2217.  
  2218.   /* return flag for lost nonzero bits */
  2219.   return ((int)j);
  2220. }
  2221. #endif
  2222.  
  2223.  
  2224. /*
  2225.  * Normalize and round off.
  2226.  *
  2227.  * The internal format number to be rounded is "s".
  2228.  * Input "lost" indicates whether or not the number is exact.
  2229.  * This is the so-called sticky bit.
  2230.  *
  2231.  * Input "subflg" indicates whether the number was obtained
  2232.  * by a subtraction operation.  In that case if lost is nonzero
  2233.  * then the number is slightly smaller than indicated.
  2234.  *
  2235.  * Input "exp" is the biased exponent, which may be negative.
  2236.  * the exponent field of "s" is ignored but is replaced by
  2237.  * "exp" as adjusted by normalization and rounding.
  2238.  *
  2239.  * Input "rcntrl" is the rounding control.
  2240.  */
  2241.  
  2242. /* For future reference:  In order for emdnorm to round off denormal
  2243.    significands at the right point, the input exponent must be
  2244.    adjusted to be the actual value it would have after conversion to
  2245.    the final floating point type.  This adjustment has been
  2246.    implemented for all type conversions (etoe53, etc.) and decimal
  2247.    conversions, but not for the arithmetic functions (eadd, etc.). 
  2248.    Data types having standard 15-bit exponents are not affected by
  2249.    this, but SFmode and DFmode are affected. For example, ediv with
  2250.    rndprc = 24 will not round correctly to 24-bit precision if the
  2251.    result is denormal.   */
  2252.  
  2253. static int rlast = -1;
  2254. static int rw = 0;
  2255. static unsigned EMUSHORT rmsk = 0;
  2256. static unsigned EMUSHORT rmbit = 0;
  2257. static unsigned EMUSHORT rebit = 0;
  2258. static int re = 0;
  2259. static unsigned EMUSHORT rbit[NI];
  2260.  
  2261. void 
  2262. emdnorm (s, lost, subflg, exp, rcntrl)
  2263.      unsigned EMUSHORT s[];
  2264.      int lost;
  2265.      int subflg;
  2266.      EMULONG exp;
  2267.      int rcntrl;
  2268. {
  2269.   int i, j;
  2270.   unsigned EMUSHORT r;
  2271.  
  2272.   /* Normalize */
  2273.   j = enormlz (s);
  2274.  
  2275.   /* a blank significand could mean either zero or infinity. */
  2276. #ifndef INFINITY
  2277.   if (j > NBITS)
  2278.     {
  2279.       ecleazs (s);
  2280.       return;
  2281.     }
  2282. #endif
  2283.   exp -= j;
  2284. #ifndef INFINITY
  2285.   if (exp >= 32767L)
  2286.     goto overf;
  2287. #else
  2288.   if ((j > NBITS) && (exp < 32767))
  2289.     {
  2290.       ecleazs (s);
  2291.       return;
  2292.     }
  2293. #endif
  2294.   if (exp < 0L)
  2295.     {
  2296.       if (exp > (EMULONG) (-NBITS - 1))
  2297.     {
  2298.       j = (int) exp;
  2299.       i = eshift (s, j);
  2300.       if (i)
  2301.         lost = 1;
  2302.     }
  2303.       else
  2304.     {
  2305.       ecleazs (s);
  2306.       return;
  2307.     }
  2308.     }
  2309.   /* Round off, unless told not to by rcntrl. */
  2310.   if (rcntrl == 0)
  2311.     goto mdfin;
  2312.   /* Set up rounding parameters if the control register changed. */
  2313.   if (rndprc != rlast)
  2314.     {
  2315.       ecleaz (rbit);
  2316.       switch (rndprc)
  2317.     {
  2318.     default:
  2319.     case NBITS:
  2320.       rw = NI - 1;        /* low guard word */
  2321.       rmsk = 0xffff;
  2322.       rmbit = 0x8000;
  2323.       re = rw - 1;
  2324.       rebit = 1;
  2325.       break;
  2326.     case 113:
  2327.       rw = 10;
  2328.       rmsk = 0x7fff;
  2329.       rmbit = 0x4000;
  2330.       rebit = 0x8000;
  2331.       re = rw;
  2332.       break;
  2333.     case 64:
  2334.       rw = 7;
  2335.       rmsk = 0xffff;
  2336.       rmbit = 0x8000;
  2337.       re = rw - 1;
  2338.       rebit = 1;
  2339.       break;
  2340.       /* For DEC or IBM arithmetic */
  2341.     case 56:
  2342.       rw = 6;
  2343.       rmsk = 0xff;
  2344.       rmbit = 0x80;
  2345.       rebit = 0x100;
  2346.       re = rw;
  2347.       break;
  2348.     case 53:
  2349.       rw = 6;
  2350.       rmsk = 0x7ff;
  2351.       rmbit = 0x0400;
  2352.       rebit = 0x800;
  2353.       re = rw;
  2354.       break;
  2355.     case 24:
  2356.       rw = 4;
  2357.       rmsk = 0xff;
  2358.       rmbit = 0x80;
  2359.       rebit = 0x100;
  2360.       re = rw;
  2361.       break;
  2362.     }
  2363.       rbit[re] = rebit;
  2364.       rlast = rndprc;
  2365.     }
  2366.  
  2367.   /* Shift down 1 temporarily if the data structure has an implied
  2368.      most significant bit and the number is denormal.  */
  2369.   if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
  2370.     {
  2371.       lost |= s[NI - 1] & 1;
  2372.       eshdn1 (s);
  2373.     }
  2374.   /* Clear out all bits below the rounding bit,
  2375.      remembering in r if any were nonzero.  */
  2376.   r = s[rw] & rmsk;
  2377.   if (rndprc < NBITS)
  2378.     {
  2379.       i = rw + 1;
  2380.       while (i < NI)
  2381.     {
  2382.       if (s[i])
  2383.         r |= 1;
  2384.       s[i] = 0;
  2385.       ++i;
  2386.     }
  2387.     }
  2388.   s[rw] &= ~rmsk;
  2389.   if ((r & rmbit) != 0)
  2390.     {
  2391.       if (r == rmbit)
  2392.     {
  2393.       if (lost == 0)
  2394.         {            /* round to even */
  2395.           if ((s[re] & rebit) == 0)
  2396.         goto mddone;
  2397.         }
  2398.       else
  2399.         {
  2400.           if (subflg != 0)
  2401.         goto mddone;
  2402.         }
  2403.     }
  2404.       eaddm (rbit, s);
  2405.     }
  2406.  mddone:
  2407.   if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
  2408.     {
  2409.       eshup1 (s);
  2410.     }
  2411.   if (s[2] != 0)
  2412.     {                /* overflow on roundoff */
  2413.       eshdn1 (s);
  2414.       exp += 1;
  2415.     }
  2416.  mdfin:
  2417.   s[NI - 1] = 0;
  2418.   if (exp >= 32767L)
  2419.     {
  2420. #ifndef INFINITY
  2421.     overf:
  2422. #endif
  2423. #ifdef INFINITY
  2424.       s[1] = 32767;
  2425.       for (i = 2; i < NI - 1; i++)
  2426.     s[i] = 0;
  2427.       if (extra_warnings)
  2428.     warning ("floating point overflow");
  2429. #else
  2430.       s[1] = 32766;
  2431.       s[2] = 0;
  2432.       for (i = M + 1; i < NI - 1; i++)
  2433.     s[i] = 0xffff;
  2434.       s[NI - 1] = 0;
  2435.       if ((rndprc < 64) || (rndprc == 113))
  2436.     {
  2437.       s[rw] &= ~rmsk;
  2438.       if (rndprc == 24)
  2439.         {
  2440.           s[5] = 0;
  2441.           s[6] = 0;
  2442.         }
  2443.     }
  2444. #endif
  2445.       return;
  2446.     }
  2447.   if (exp < 0)
  2448.     s[1] = 0;
  2449.   else
  2450.     s[1] = (unsigned EMUSHORT) exp;
  2451. }
  2452.  
  2453.  
  2454.  
  2455. /*
  2456. ;    Subtract external format numbers.
  2457. ;
  2458. ;    unsigned EMUSHORT a[NE], b[NE], c[NE];
  2459. ;    esub (a, b, c);     c = b - a
  2460. */
  2461.  
  2462. static int subflg = 0;
  2463.  
  2464. void 
  2465. esub (a, b, c)
  2466.      unsigned EMUSHORT *a, *b, *c;
  2467. {
  2468.  
  2469. #ifdef NANS
  2470.   if (eisnan (a))
  2471.     {
  2472.       emov (a, c);
  2473.       return;
  2474.     }
  2475.   if (eisnan (b))
  2476.     {
  2477.       emov (b, c);
  2478.       return;
  2479.     }
  2480. /* Infinity minus infinity is a NaN.
  2481.    Test for subtracting infinities of the same sign. */
  2482.   if (eisinf (a) && eisinf (b)
  2483.       && ((eisneg (a) ^ eisneg (b)) == 0))
  2484.     {
  2485.       mtherr ("esub", INVALID);
  2486.       enan (c);
  2487.       return;
  2488.     }
  2489. #endif
  2490.   subflg = 1;
  2491.   eadd1 (a, b, c);
  2492. }
  2493.  
  2494.  
  2495. /*
  2496. ;    Add.
  2497. ;
  2498. ;    unsigned EMUSHORT a[NE], b[NE], c[NE];
  2499. ;    eadd (a, b, c);     c = b + a
  2500. */
  2501. void 
  2502. eadd (a, b, c)
  2503.      unsigned EMUSHORT *a, *b, *c;
  2504. {
  2505.  
  2506. #ifdef NANS
  2507. /* NaN plus anything is a NaN. */
  2508.   if (eisnan (a))
  2509.     {
  2510.       emov (a, c);
  2511.       return;
  2512.     }
  2513.   if (eisnan (b))
  2514.     {
  2515.       emov (b, c);
  2516.       return;
  2517.     }
  2518. /* Infinity minus infinity is a NaN.
  2519.    Test for adding infinities of opposite signs. */
  2520.   if (eisinf (a) && eisinf (b)
  2521.       && ((eisneg (a) ^ eisneg (b)) != 0))
  2522.     {
  2523.       mtherr ("esub", INVALID);
  2524.       enan (c);
  2525.       return;
  2526.     }
  2527. #endif
  2528.   subflg = 0;
  2529.   eadd1 (a, b, c);
  2530. }
  2531.  
  2532. void 
  2533. eadd1 (a, b, c)
  2534.      unsigned EMUSHORT *a, *b, *c;
  2535. {
  2536.   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
  2537.   int i, lost, j, k;
  2538.   EMULONG lt, lta, ltb;
  2539.  
  2540. #ifdef INFINITY
  2541.   if (eisinf (a))
  2542.     {
  2543.       emov (a, c);
  2544.       if (subflg)
  2545.     eneg (c);
  2546.       return;
  2547.     }
  2548.   if (eisinf (b))
  2549.     {
  2550.       emov (b, c);
  2551.       return;
  2552.     }
  2553. #endif
  2554.   emovi (a, ai);
  2555.   emovi (b, bi);
  2556.   if (subflg)
  2557.     ai[0] = ~ai[0];
  2558.  
  2559.   /* compare exponents */
  2560.   lta = ai[E];
  2561.   ltb = bi[E];
  2562.   lt = lta - ltb;
  2563.   if (lt > 0L)
  2564.     {                /* put the larger number in bi */
  2565.       emovz (bi, ci);
  2566.       emovz (ai, bi);
  2567.       emovz (ci, ai);
  2568.       ltb = bi[E];
  2569.       lt = -lt;
  2570.     }
  2571.   lost = 0;
  2572.   if (lt != 0L)
  2573.     {
  2574.       if (lt < (EMULONG) (-NBITS - 1))
  2575.     goto done;        /* answer same as larger addend */
  2576.       k = (int) lt;
  2577.       lost = eshift (ai, k);    /* shift the smaller number down */
  2578.     }
  2579.   else
  2580.     {
  2581.       /* exponents were the same, so must compare significands */
  2582.       i = ecmpm (ai, bi);
  2583.       if (i == 0)
  2584.     {            /* the numbers are identical in magnitude */
  2585.       /* if different signs, result is zero */
  2586.       if (ai[0] != bi[0])
  2587.         {
  2588.           eclear (c);
  2589.           return;
  2590.         }
  2591.       /* if same sign, result is double */
  2592.       /* double denomalized tiny number */
  2593.       if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
  2594.         {
  2595.           eshup1 (bi);
  2596.           goto done;
  2597.         }
  2598.       /* add 1 to exponent unless both are zero! */
  2599.       for (j = 1; j < NI - 1; j++)
  2600.         {
  2601.           if (bi[j] != 0)
  2602.         {
  2603.           /* This could overflow, but let emovo take care of that. */
  2604.           ltb += 1;
  2605.           break;
  2606.         }
  2607.         }
  2608.       bi[E] = (unsigned EMUSHORT) ltb;
  2609.       goto done;
  2610.     }
  2611.       if (i > 0)
  2612.     {            /* put the larger number in bi */
  2613.       emovz (bi, ci);
  2614.       emovz (ai, bi);
  2615.       emovz (ci, ai);
  2616.     }
  2617.     }
  2618.   if (ai[0] == bi[0])
  2619.     {
  2620.       eaddm (ai, bi);
  2621.       subflg = 0;
  2622.     }
  2623.   else
  2624.     {
  2625.       esubm (ai, bi);
  2626.       subflg = 1;
  2627.     }
  2628.   emdnorm (bi, lost, subflg, ltb, 64);
  2629.  
  2630.  done:
  2631.   emovo (bi, c);
  2632. }
  2633.  
  2634.  
  2635.  
  2636. /*
  2637. ;    Divide.
  2638. ;
  2639. ;    unsigned EMUSHORT a[NE], b[NE], c[NE];
  2640. ;    ediv (a, b, c);    c = b / a
  2641. */
  2642. void 
  2643. ediv (a, b, c)
  2644.      unsigned EMUSHORT *a, *b, *c;
  2645. {
  2646.   unsigned EMUSHORT ai[NI], bi[NI];
  2647.   int i;
  2648.   EMULONG lt, lta, ltb;
  2649.  
  2650. #ifdef NANS
  2651. /* Return any NaN input. */
  2652.   if (eisnan (a))
  2653.     {
  2654.     emov (a, c);
  2655.     return;
  2656.     }
  2657.   if (eisnan (b))
  2658.     {
  2659.     emov (b, c);
  2660.     return;
  2661.     }
  2662. /* Zero over zero, or infinity over infinity, is a NaN. */
  2663.   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
  2664.       || (eisinf (a) && eisinf (b)))
  2665.     {
  2666.     mtherr ("ediv", INVALID);
  2667.     enan (c);
  2668.     return;
  2669.     }
  2670. #endif
  2671. /* Infinity over anything else is infinity. */
  2672. #ifdef INFINITY
  2673.   if (eisinf (b))
  2674.     {
  2675.       if (eisneg (a) ^ eisneg (b))
  2676.     *(c + (NE - 1)) = 0x8000;
  2677.       else
  2678.     *(c + (NE - 1)) = 0;
  2679.       einfin (c);
  2680.       return;
  2681.     }
  2682. /* Anything else over infinity is zero. */
  2683.   if (eisinf (a))
  2684.     {
  2685.       eclear (c);
  2686.       return;
  2687.     }
  2688. #endif
  2689.   emovi (a, ai);
  2690.   emovi (b, bi);
  2691.   lta = ai[E];
  2692.   ltb = bi[E];
  2693.   if (bi[E] == 0)
  2694.     {                /* See if numerator is zero. */
  2695.       for (i = 1; i < NI - 1; i++)
  2696.     {
  2697.       if (bi[i] != 0)
  2698.         {
  2699.           ltb -= enormlz (bi);
  2700.           goto dnzro1;
  2701.         }
  2702.     }
  2703.       eclear (c);
  2704.       return;
  2705.     }
  2706.  dnzro1:
  2707.  
  2708.   if (ai[E] == 0)
  2709.     {                /* possible divide by zero */
  2710.       for (i = 1; i < NI - 1; i++)
  2711.     {
  2712.       if (ai[i] != 0)
  2713.         {
  2714.           lta -= enormlz (ai);
  2715.           goto dnzro2;
  2716.         }
  2717.     }
  2718.       if (ai[0] == bi[0])
  2719.     *(c + (NE - 1)) = 0;
  2720.       else
  2721.     *(c + (NE - 1)) = 0x8000;
  2722. /* Divide by zero is not an invalid operation.
  2723.    It is a divide-by-zero operation!   */
  2724.       einfin (c);
  2725.       mtherr ("ediv", SING);
  2726.       return;
  2727.     }
  2728.  dnzro2:
  2729.  
  2730.   i = edivm (ai, bi);
  2731.   /* calculate exponent */
  2732.   lt = ltb - lta + EXONE;
  2733.   emdnorm (bi, i, 0, lt, 64);
  2734.   /* set the sign */
  2735.   if (ai[0] == bi[0])
  2736.     bi[0] = 0;
  2737.   else
  2738.     bi[0] = 0Xffff;
  2739.   emovo (bi, c);
  2740. }
  2741.  
  2742.  
  2743.  
  2744. /*
  2745. ;    Multiply.
  2746. ;
  2747. ;    unsigned EMUSHORT a[NE], b[NE], c[NE];
  2748. ;    emul (a, b, c);    c = b * a
  2749. */
  2750. void 
  2751. emul (a, b, c)
  2752.      unsigned EMUSHORT *a, *b, *c;
  2753. {
  2754.   unsigned EMUSHORT ai[NI], bi[NI];
  2755.   int i, j;
  2756.   EMULONG lt, lta, ltb;
  2757.  
  2758. #ifdef NANS
  2759. /* NaN times anything is the same NaN. */
  2760.   if (eisnan (a))
  2761.     {
  2762.     emov (a, c);
  2763.     return;
  2764.     }
  2765.   if (eisnan (b))
  2766.     {
  2767.     emov (b, c);
  2768.     return;
  2769.     }
  2770. /* Zero times infinity is a NaN. */
  2771.   if ((eisinf (a) && (ecmp (b, ezero) == 0))
  2772.       || (eisinf (b) && (ecmp (a, ezero) == 0)))
  2773.     {
  2774.     mtherr ("emul", INVALID);
  2775.     enan (c);
  2776.     return;
  2777.     }
  2778. #endif
  2779. /* Infinity times anything else is infinity. */
  2780. #ifdef INFINITY
  2781.   if (eisinf (a) || eisinf (b))
  2782.     {
  2783.       if (eisneg (a) ^ eisneg (b))
  2784.     *(c + (NE - 1)) = 0x8000;
  2785.       else
  2786.     *(c + (NE - 1)) = 0;
  2787.       einfin (c);
  2788.       return;
  2789.     }
  2790. #endif
  2791.   emovi (a, ai);
  2792.   emovi (b, bi);
  2793.   lta = ai[E];
  2794.   ltb = bi[E];
  2795.   if (ai[E] == 0)
  2796.     {
  2797.       for (i = 1; i < NI - 1; i++)
  2798.     {
  2799.       if (ai[i] != 0)
  2800.         {
  2801.           lta -= enormlz (ai);
  2802.           goto mnzer1;
  2803.         }
  2804.     }
  2805.       eclear (c);
  2806.       return;
  2807.     }
  2808.  mnzer1:
  2809.  
  2810.   if (bi[E] == 0)
  2811.     {
  2812.       for (i = 1; i < NI - 1; i++)
  2813.     {
  2814.       if (bi[i] != 0)
  2815.         {
  2816.           ltb -= enormlz (bi);
  2817.           goto mnzer2;
  2818.         }
  2819.     }
  2820.       eclear (c);
  2821.       return;
  2822.     }
  2823.  mnzer2:
  2824.  
  2825.   /* Multiply significands */
  2826.   j = emulm (ai, bi);
  2827.   /* calculate exponent */
  2828.   lt = lta + ltb - (EXONE - 1);
  2829.   emdnorm (bi, j, 0, lt, 64);
  2830.   /* calculate sign of product */
  2831.   if (ai[0] == bi[0])
  2832.     bi[0] = 0;
  2833.   else
  2834.     bi[0] = 0xffff;
  2835.   emovo (bi, c);
  2836. }
  2837.  
  2838.  
  2839.  
  2840.  
  2841. /*
  2842. ; Convert IEEE double precision to e type
  2843. ;    double d;
  2844. ;    unsigned EMUSHORT x[N+2];
  2845. ;    e53toe (&d, x);
  2846. */
  2847. void 
  2848. e53toe (pe, y)
  2849.      unsigned EMUSHORT *pe, *y;
  2850. {
  2851. #ifdef DEC
  2852.  
  2853.   dectoe (pe, y);        /* see etodec.c */
  2854.  
  2855. #else
  2856. #ifdef IBM
  2857.  
  2858.   ibmtoe (pe, y, DFmode);
  2859.  
  2860. #else
  2861.   register unsigned EMUSHORT r;
  2862.   register unsigned EMUSHORT *e, *p;
  2863.   unsigned EMUSHORT yy[NI];
  2864.   int denorm, k;
  2865.  
  2866.   e = pe;
  2867.   denorm = 0;            /* flag if denormalized number */
  2868.   ecleaz (yy);
  2869. #ifdef IBMPC
  2870.   e += 3;
  2871. #endif
  2872.   r = *e;
  2873.   yy[0] = 0;
  2874.   if (r & 0x8000)
  2875.     yy[0] = 0xffff;
  2876.   yy[M] = (r & 0x0f) | 0x10;
  2877.   r &= ~0x800f;            /* strip sign and 4 significand bits */
  2878. #ifdef INFINITY
  2879.   if (r == 0x7ff0)
  2880.     {
  2881. #ifdef NANS
  2882. #ifdef IBMPC
  2883.       if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
  2884.       || (pe[1] != 0) || (pe[0] != 0))
  2885.     {
  2886.       enan (y);
  2887.       return;
  2888.     }
  2889. #else
  2890.       if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
  2891.       || (pe[2] != 0) || (pe[3] != 0))
  2892.     {
  2893.       enan (y);
  2894.       return;
  2895.     }
  2896. #endif
  2897. #endif  /* NANS */
  2898.       eclear (y);
  2899.       einfin (y);
  2900.       if (yy[0])
  2901.     eneg (y);
  2902.       return;
  2903.     }
  2904. #endif  /* INFINITY */
  2905.   r >>= 4;
  2906.   /* If zero exponent, then the significand is denormalized.
  2907.    * So, take back the understood high significand bit. */
  2908.   if (r == 0)
  2909.     {
  2910.       denorm = 1;
  2911.       yy[M] &= ~0x10;
  2912.     }
  2913.   r += EXONE - 01777;
  2914.   yy[E] = r;
  2915.   p = &yy[M + 1];
  2916. #ifdef IBMPC
  2917.   *p++ = *(--e);
  2918.   *p++ = *(--e);
  2919.   *p++ = *(--e);
  2920. #endif
  2921. #ifdef MIEEE
  2922.   ++e;
  2923.   *p++ = *e++;
  2924.   *p++ = *e++;
  2925.   *p++ = *e++;
  2926. #endif
  2927.   eshift (yy, -5);
  2928.   if (denorm)
  2929.     {                /* if zero exponent, then normalize the significand */
  2930.       if ((k = enormlz (yy)) > NBITS)
  2931.     ecleazs (yy);
  2932.       else
  2933.     yy[E] -= (unsigned EMUSHORT) (k - 1);
  2934.     }
  2935.   emovo (yy, y);
  2936. #endif /* not IBM */
  2937. #endif /* not DEC */
  2938. }
  2939.  
  2940. void 
  2941. e64toe (pe, y)
  2942.      unsigned EMUSHORT *pe, *y;
  2943. {
  2944.   unsigned EMUSHORT yy[NI];
  2945.   unsigned EMUSHORT *e, *p, *q;
  2946.   int i;
  2947.  
  2948.   e = pe;
  2949.   p = yy;
  2950.   for (i = 0; i < NE - 5; i++)
  2951.     *p++ = 0;
  2952. #ifdef IBMPC
  2953.   for (i = 0; i < 5; i++)
  2954.     *p++ = *e++;
  2955. #endif
  2956. /* This precision is not ordinarily supported on DEC or IBM. */
  2957. #ifdef DEC
  2958.   for (i = 0; i < 5; i++)
  2959.     *p++ = *e++;
  2960. #endif
  2961. #ifdef IBM
  2962.   p = &yy[0] + (NE - 1);
  2963.   *p-- = *e++;
  2964.   ++e;
  2965.   for (i = 0; i < 5; i++)
  2966.     *p-- = *e++;
  2967. #endif
  2968. #ifdef MIEEE
  2969.   p = &yy[0] + (NE - 1);
  2970.   *p-- = *e++;
  2971.   ++e;
  2972.   for (i = 0; i < 4; i++)
  2973.     *p-- = *e++;
  2974. #endif
  2975.   p = yy;
  2976.   q = y;
  2977. #ifdef INFINITY
  2978.   if (*p == 0x7fff)
  2979.     {
  2980. #ifdef NANS
  2981. #ifdef IBMPC
  2982.       for (i = 0; i < 4; i++)
  2983.     {
  2984.       if (pe[i] != 0)
  2985.         {
  2986.           enan (y);
  2987.           return;
  2988.         }
  2989.     }
  2990. #else
  2991.       for (i = 1; i <= 4; i++)
  2992.     {
  2993.       if (pe[i] != 0)
  2994.         {
  2995.           enan (y);
  2996.           return;
  2997.         }
  2998.     }
  2999. #endif
  3000. #endif /* NANS */
  3001.       eclear (y);
  3002.       einfin (y);
  3003.       if (*p & 0x8000)
  3004.     eneg (y);
  3005.       return;
  3006.     }
  3007. #endif  /* INFINITY */
  3008.   for (i = 0; i < NE; i++)
  3009.     *q++ = *p++;
  3010. }
  3011.  
  3012.  
  3013. void 
  3014. e113toe (pe, y)
  3015.      unsigned EMUSHORT *pe, *y;
  3016. {
  3017.   register unsigned EMUSHORT r;
  3018.   unsigned EMUSHORT *e, *p;
  3019.   unsigned EMUSHORT yy[NI];
  3020.   int denorm, i;
  3021.  
  3022.   e = pe;
  3023.   denorm = 0;
  3024.   ecleaz (yy);
  3025. #ifdef IBMPC
  3026.   e += 7;
  3027. #endif
  3028.   r = *e;
  3029.   yy[0] = 0;
  3030.   if (r & 0x8000)
  3031.     yy[0] = 0xffff;
  3032.   r &= 0x7fff;
  3033. #ifdef INFINITY
  3034.   if (r == 0x7fff)
  3035.     {
  3036. #ifdef NANS
  3037. #ifdef IBMPC
  3038.       for (i = 0; i < 7; i++)
  3039.     {
  3040.       if (pe[i] != 0)
  3041.         {
  3042.           enan (y);
  3043.           return;
  3044.         }
  3045.     }
  3046. #else
  3047.       for (i = 1; i < 8; i++)
  3048.     {
  3049.       if (pe[i] != 0)
  3050.         {
  3051.           enan (y);
  3052.           return;
  3053.         }
  3054.     }
  3055. #endif
  3056. #endif /* NANS */
  3057.       eclear (y);
  3058.       einfin (y);
  3059.       if (yy[0])
  3060.     eneg (y);
  3061.       return;
  3062.     }
  3063. #endif  /* INFINITY */
  3064.   yy[E] = r;
  3065.   p = &yy[M + 1];
  3066. #ifdef IBMPC
  3067.   for (i = 0; i < 7; i++)
  3068.     *p++ = *(--e);
  3069. #endif
  3070. #ifdef MIEEE
  3071.   ++e;
  3072.   for (i = 0; i < 7; i++)
  3073.     *p++ = *e++;
  3074. #endif
  3075. /* If denormal, remove the implied bit; else shift down 1. */
  3076.   if (r == 0)
  3077.     {
  3078.       yy[M] = 0;
  3079.     }
  3080.   else
  3081.     {
  3082.       yy[M] = 1;
  3083.       eshift (yy, -1);
  3084.     }
  3085.   emovo (yy, y);
  3086. }
  3087.  
  3088.  
  3089. /*
  3090. ; Convert IEEE single precision to e type
  3091. ;    float d;
  3092. ;    unsigned EMUSHORT x[N+2];
  3093. ;    dtox (&d, x);
  3094. */
  3095. void 
  3096. e24toe (pe, y)
  3097.      unsigned EMUSHORT *pe, *y;
  3098. {
  3099. #ifdef IBM
  3100.  
  3101.   ibmtoe (pe, y, SFmode);
  3102.  
  3103. #else
  3104.   register unsigned EMUSHORT r;
  3105.   register unsigned EMUSHORT *e, *p;
  3106.   unsigned EMUSHORT yy[NI];
  3107.   int denorm, k;
  3108.  
  3109.   e = pe;
  3110.   denorm = 0;            /* flag if denormalized number */
  3111.   ecleaz (yy);
  3112. #ifdef IBMPC
  3113.   e += 1;
  3114. #endif
  3115. #ifdef DEC
  3116.   e += 1;
  3117. #endif
  3118.   r = *e;
  3119.   yy[0] = 0;
  3120.   if (r & 0x8000)
  3121.     yy[0] = 0xffff;
  3122.   yy[M] = (r & 0x7f) | 0200;
  3123.   r &= ~0x807f;            /* strip sign and 7 significand bits */
  3124. #ifdef INFINITY
  3125.   if (r == 0x7f80)
  3126.     {
  3127. #ifdef NANS
  3128. #ifdef MIEEE
  3129.       if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
  3130.     {
  3131.       enan (y);
  3132.       return;
  3133.     }
  3134. #else
  3135.       if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
  3136.     {
  3137.       enan (y);
  3138.       return;
  3139.     }
  3140. #endif
  3141. #endif  /* NANS */
  3142.       eclear (y);
  3143.       einfin (y);
  3144.       if (yy[0])
  3145.     eneg (y);
  3146.       return;
  3147.     }
  3148. #endif  /* INFINITY */
  3149.   r >>= 7;
  3150.   /* If zero exponent, then the significand is denormalized.
  3151.    * So, take back the understood high significand bit. */
  3152.   if (r == 0)
  3153.     {
  3154.       denorm = 1;
  3155.       yy[M] &= ~0200;
  3156.     }
  3157.   r += EXONE - 0177;
  3158.   yy[E] = r;
  3159.   p = &yy[M + 1];
  3160. #ifdef IBMPC
  3161.   *p++ = *(--e);
  3162. #endif
  3163. #ifdef DEC
  3164.   *p++ = *(--e);
  3165. #endif
  3166. #ifdef MIEEE
  3167.   ++e;
  3168.   *p++ = *e++;
  3169. #endif
  3170.   eshift (yy, -8);
  3171.   if (denorm)
  3172.     {                /* if zero exponent, then normalize the significand */
  3173.       if ((k = enormlz (yy)) > NBITS)
  3174.     ecleazs (yy);
  3175.       else
  3176.     yy[E] -= (unsigned EMUSHORT) (k - 1);
  3177.     }
  3178.   emovo (yy, y);
  3179. #endif /* not IBM */
  3180. }
  3181.  
  3182.  
  3183. void 
  3184. etoe113 (x, e)
  3185.      unsigned EMUSHORT *x, *e;
  3186. {
  3187.   unsigned EMUSHORT xi[NI];
  3188.   EMULONG exp;
  3189.   int rndsav;
  3190.  
  3191. #ifdef NANS
  3192.   if (eisnan (x))
  3193.     {
  3194.       make_nan (e, TFmode);
  3195.       return;
  3196.     }
  3197. #endif
  3198.   emovi (x, xi);
  3199.   exp = (EMULONG) xi[E];
  3200. #ifdef INFINITY
  3201.   if (eisinf (x))
  3202.     goto nonorm;
  3203. #endif
  3204.   /* round off to nearest or even */
  3205.   rndsav = rndprc;
  3206.   rndprc = 113;
  3207.   emdnorm (xi, 0, 0, exp, 64);
  3208.   rndprc = rndsav;
  3209.  nonorm:
  3210.   toe113 (xi, e);
  3211. }
  3212.  
  3213. /* move out internal format to ieee long double */
  3214. static void 
  3215. toe113 (a, b)
  3216.      unsigned EMUSHORT *a, *b;
  3217. {
  3218.   register unsigned EMUSHORT *p, *q;
  3219.   unsigned EMUSHORT i;
  3220.  
  3221. #ifdef NANS
  3222.   if (eiisnan (a))
  3223.     {
  3224.       make_nan (b, TFmode);
  3225.       return;
  3226.     }
  3227. #endif
  3228.   p = a;
  3229. #ifdef MIEEE
  3230.   q = b;
  3231. #else
  3232.   q = b + 7;            /* point to output exponent */
  3233. #endif
  3234.  
  3235.   /* If not denormal, delete the implied bit. */
  3236.   if (a[E] != 0)
  3237.     {
  3238.       eshup1 (a);
  3239.     }
  3240.   /* combine sign and exponent */
  3241.   i = *p++;
  3242. #ifdef MIEEE
  3243.   if (i)
  3244.     *q++ = *p++ | 0x8000;
  3245.   else
  3246.     *q++ = *p++;
  3247. #else
  3248.   if (i)
  3249.     *q-- = *p++ | 0x8000;
  3250.   else
  3251.     *q-- = *p++;
  3252. #endif
  3253.   /* skip over guard word */
  3254.   ++p;
  3255.   /* move the significand */
  3256. #ifdef MIEEE
  3257.   for (i = 0; i < 7; i++)
  3258.     *q++ = *p++;
  3259. #else
  3260.   for (i = 0; i < 7; i++)
  3261.     *q-- = *p++;
  3262. #endif
  3263. }
  3264.  
  3265. void 
  3266. etoe64 (x, e)
  3267.      unsigned EMUSHORT *x, *e;
  3268. {
  3269.   unsigned EMUSHORT xi[NI];
  3270.   EMULONG exp;
  3271.   int rndsav;
  3272.  
  3273. #ifdef NANS
  3274.   if (eisnan (x))
  3275.     {
  3276.       make_nan (e, XFmode);
  3277.       return;
  3278.     }
  3279. #endif
  3280.   emovi (x, xi);
  3281.   /* adjust exponent for offset */
  3282.   exp = (EMULONG) xi[E];
  3283. #ifdef INFINITY
  3284.   if (eisinf (x))
  3285.     goto nonorm;
  3286. #endif
  3287.   /* round off to nearest or even */
  3288.   rndsav = rndprc;
  3289.   rndprc = 64;
  3290.   emdnorm (xi, 0, 0, exp, 64);
  3291.   rndprc = rndsav;
  3292.  nonorm:
  3293.   toe64 (xi, e);
  3294. }
  3295.  
  3296. /* move out internal format to ieee long double */
  3297. static void 
  3298. toe64 (a, b)
  3299.      unsigned EMUSHORT *a, *b;
  3300. {
  3301.   register unsigned EMUSHORT *p, *q;
  3302.   unsigned EMUSHORT i;
  3303.  
  3304. #ifdef NANS
  3305.   if (eiisnan (a))
  3306.     {
  3307.       make_nan (b, XFmode);
  3308.       return;
  3309.     }
  3310. #endif
  3311.   p = a;
  3312. #if defined(MIEEE) || defined(IBM)
  3313.   q = b;
  3314. #else
  3315.   q = b + 4;            /* point to output exponent */
  3316. #if LONG_DOUBLE_TYPE_SIZE == 96
  3317.   /* Clear the last two bytes of 12-byte Intel format */
  3318.   *(q+1) = 0;
  3319. #endif
  3320. #endif
  3321.  
  3322.   /* combine sign and exponent */
  3323.   i = *p++;
  3324. #if defined(MIEEE) || defined(IBM)
  3325.   if (i)
  3326.     *q++ = *p++ | 0x8000;
  3327.   else
  3328.     *q++ = *p++;
  3329.   *q++ = 0;
  3330. #else
  3331.   if (i)
  3332.     *q-- = *p++ | 0x8000;
  3333.   else
  3334.     *q-- = *p++;
  3335. #endif
  3336.   /* skip over guard word */
  3337.   ++p;
  3338.   /* move the significand */
  3339. #if defined(MIEEE) || defined(IBM)
  3340.   for (i = 0; i < 4; i++)
  3341.     *q++ = *p++;
  3342. #else
  3343.   for (i = 0; i < 4; i++)
  3344.     *q-- = *p++;
  3345. #endif
  3346. }
  3347.  
  3348.  
  3349. /*
  3350. ; e type to IEEE double precision
  3351. ;    double d;
  3352. ;    unsigned EMUSHORT x[NE];
  3353. ;    etoe53 (x, &d);
  3354. */
  3355.  
  3356. #ifdef DEC
  3357.  
  3358. void 
  3359. etoe53 (x, e)
  3360.      unsigned EMUSHORT *x, *e;
  3361. {
  3362.   etodec (x, e);        /* see etodec.c */
  3363. }
  3364.  
  3365. static void 
  3366. toe53 (x, y)
  3367.      unsigned EMUSHORT *x, *y;
  3368. {
  3369.   todec (x, y);
  3370. }
  3371.  
  3372. #else
  3373. #ifdef IBM
  3374.  
  3375. void 
  3376. etoe53 (x, e)
  3377.      unsigned EMUSHORT *x, *e;
  3378. {
  3379.   etoibm (x, e, DFmode);
  3380. }
  3381.  
  3382. static void 
  3383. toe53 (x, y)
  3384.      unsigned EMUSHORT *x, *y;
  3385. {
  3386.   toibm (x, y, DFmode);
  3387. }
  3388.  
  3389. #else  /* it's neither DEC nor IBM */
  3390.  
  3391. void 
  3392. etoe53 (x, e)
  3393.      unsigned EMUSHORT *x, *e;
  3394. {
  3395.   unsigned EMUSHORT xi[NI];
  3396.   EMULONG exp;
  3397.   int rndsav;
  3398.  
  3399. #ifdef NANS
  3400.   if (eisnan (x))
  3401.     {
  3402.       make_nan (e, DFmode);
  3403.       return;
  3404.     }
  3405. #endif
  3406.   emovi (x, xi);
  3407.   /* adjust exponent for offsets */
  3408.   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
  3409. #ifdef INFINITY
  3410.   if (eisinf (x))
  3411.     goto nonorm;
  3412. #endif
  3413.   /* round off to nearest or even */
  3414.   rndsav = rndprc;
  3415.   rndprc = 53;
  3416.   emdnorm (xi, 0, 0, exp, 64);
  3417.   rndprc = rndsav;
  3418.  nonorm:
  3419.   toe53 (xi, e);
  3420. }
  3421.  
  3422.  
  3423. static void 
  3424. toe53 (x, y)
  3425.      unsigned EMUSHORT *x, *y;
  3426. {
  3427.   unsigned EMUSHORT i;
  3428.   unsigned EMUSHORT *p;
  3429.  
  3430. #ifdef NANS
  3431.   if (eiisnan (x))
  3432.     {
  3433.       make_nan (y, DFmode);
  3434.       return;
  3435.     }
  3436. #endif
  3437.   p = &x[0];
  3438. #ifdef IBMPC
  3439.   y += 3;
  3440. #endif
  3441.   *y = 0;            /* output high order */
  3442.   if (*p++)
  3443.     *y = 0x8000;        /* output sign bit */
  3444.  
  3445.   i = *p++;
  3446.   if (i >= (unsigned int) 2047)
  3447.     {                /* Saturate at largest number less than infinity. */
  3448. #ifdef INFINITY
  3449.       *y |= 0x7ff0;
  3450. #ifdef IBMPC
  3451.       *(--y) = 0;
  3452.       *(--y) = 0;
  3453.       *(--y) = 0;
  3454. #endif
  3455. #ifdef MIEEE
  3456.       ++y;
  3457.       *y++ = 0;
  3458.       *y++ = 0;
  3459.       *y++ = 0;
  3460. #endif
  3461. #else
  3462.       *y |= (unsigned EMUSHORT) 0x7fef;
  3463. #ifdef IBMPC
  3464.       *(--y) = 0xffff;
  3465.       *(--y) = 0xffff;
  3466.       *(--y) = 0xffff;
  3467. #endif
  3468. #ifdef MIEEE
  3469.       ++y;
  3470.       *y++ = 0xffff;
  3471.       *y++ = 0xffff;
  3472.       *y++ = 0xffff;
  3473. #endif
  3474. #endif
  3475.       return;
  3476.     }
  3477.   if (i == 0)
  3478.     {
  3479.       eshift (x, 4);
  3480.     }
  3481.   else
  3482.     {
  3483.       i <<= 4;
  3484.       eshift (x, 5);
  3485.     }
  3486.   i |= *p++ & (unsigned EMUSHORT) 0x0f;    /* *p = xi[M] */
  3487.   *y |= (unsigned EMUSHORT) i;    /* high order output already has sign bit set */
  3488. #ifdef IBMPC
  3489.   *(--y) = *p++;
  3490.   *(--y) = *p++;
  3491.   *(--y) = *p;
  3492. #endif
  3493. #ifdef MIEEE
  3494.   ++y;
  3495.   *y++ = *p++;
  3496.   *y++ = *p++;
  3497.   *y++ = *p++;
  3498. #endif
  3499. }
  3500.  
  3501. #endif /* not IBM */
  3502. #endif /* not DEC */
  3503.  
  3504.  
  3505.  
  3506. /*
  3507. ; e type to IEEE single precision
  3508. ;    float d;
  3509. ;    unsigned EMUSHORT x[N+2];
  3510. ;    xtod (x, &d);
  3511. */
  3512. #ifdef IBM
  3513.  
  3514. void 
  3515. etoe24 (x, e)
  3516.      unsigned EMUSHORT *x, *e;
  3517. {
  3518.   etoibm (x, e, SFmode);
  3519. }
  3520.  
  3521. static void 
  3522. toe24 (x, y)
  3523.      unsigned EMUSHORT *x, *y;
  3524. {
  3525.   toibm (x, y, SFmode);
  3526. }
  3527.  
  3528. #else
  3529.  
  3530. void 
  3531. etoe24 (x, e)
  3532.      unsigned EMUSHORT *x, *e;
  3533. {
  3534.   EMULONG exp;
  3535.   unsigned EMUSHORT xi[NI];
  3536.   int rndsav;
  3537.  
  3538. #ifdef NANS
  3539.   if (eisnan (x))
  3540.     {
  3541.       make_nan (e, SFmode);
  3542.       return;
  3543.     }
  3544. #endif
  3545.   emovi (x, xi);
  3546.   /* adjust exponent for offsets */
  3547.   exp = (EMULONG) xi[E] - (EXONE - 0177);
  3548. #ifdef INFINITY
  3549.   if (eisinf (x))
  3550.     goto nonorm;
  3551. #endif
  3552.   /* round off to nearest or even */
  3553.   rndsav = rndprc;
  3554.   rndprc = 24;
  3555.   emdnorm (xi, 0, 0, exp, 64);
  3556.   rndprc = rndsav;
  3557.  nonorm:
  3558.   toe24 (xi, e);
  3559. }
  3560.  
  3561. static void 
  3562. toe24 (x, y)
  3563.      unsigned EMUSHORT *x, *y;
  3564. {
  3565.   unsigned EMUSHORT i;
  3566.   unsigned EMUSHORT *p;
  3567.  
  3568. #ifdef NANS
  3569.   if (eiisnan (x))
  3570.     {
  3571.       make_nan (y, SFmode);
  3572.       return;
  3573.     }
  3574. #endif
  3575.   p = &x[0];
  3576. #ifdef IBMPC
  3577.   y += 1;
  3578. #endif
  3579. #ifdef DEC
  3580.   y += 1;
  3581. #endif
  3582.   *y = 0;            /* output high order */
  3583.   if (*p++)
  3584.     *y = 0x8000;        /* output sign bit */
  3585.  
  3586.   i = *p++;
  3587. /* Handle overflow cases. */
  3588.   if (i >= 255)
  3589.     {
  3590. #ifdef INFINITY
  3591.       *y |= (unsigned EMUSHORT) 0x7f80;
  3592. #ifdef IBMPC
  3593.       *(--y) = 0;
  3594. #endif
  3595. #ifdef DEC
  3596.       *(--y) = 0;
  3597. #endif
  3598. #ifdef MIEEE
  3599.       ++y;
  3600.       *y = 0;
  3601. #endif
  3602. #else  /* no INFINITY */
  3603.       *y |= (unsigned EMUSHORT) 0x7f7f;
  3604. #ifdef IBMPC
  3605.       *(--y) = 0xffff;
  3606. #endif
  3607. #ifdef DEC
  3608.       *(--y) = 0xffff;
  3609. #endif
  3610. #ifdef MIEEE
  3611.       ++y;
  3612.       *y = 0xffff;
  3613. #endif
  3614. #ifdef ERANGE
  3615.       errno = ERANGE;
  3616. #endif
  3617. #endif  /* no INFINITY */
  3618.       return;
  3619.     }
  3620.   if (i == 0)
  3621.     {
  3622.       eshift (x, 7);
  3623.     }
  3624.   else
  3625.     {
  3626.       i <<= 7;
  3627.       eshift (x, 8);
  3628.     }
  3629.   i |= *p++ & (unsigned EMUSHORT) 0x7f;    /* *p = xi[M] */
  3630.   *y |= i;            /* high order output already has sign bit set */
  3631. #ifdef IBMPC
  3632.   *(--y) = *p;
  3633. #endif
  3634. #ifdef DEC
  3635.   *(--y) = *p;
  3636. #endif
  3637. #ifdef MIEEE
  3638.   ++y;
  3639.   *y = *p;
  3640. #endif
  3641. }
  3642. #endif  /* not IBM */
  3643.  
  3644. /* Compare two e type numbers.
  3645.  *
  3646.  * unsigned EMUSHORT a[NE], b[NE];
  3647.  * ecmp (a, b);
  3648.  *
  3649.  *  returns +1 if a > b
  3650.  *           0 if a == b
  3651.  *          -1 if a < b
  3652.  *          -2 if either a or b is a NaN.
  3653.  */
  3654. int 
  3655. ecmp (a, b)
  3656.      unsigned EMUSHORT *a, *b;
  3657. {
  3658.   unsigned EMUSHORT ai[NI], bi[NI];
  3659.   register unsigned EMUSHORT *p, *q;
  3660.   register int i;
  3661.   int msign;
  3662.  
  3663. #ifdef NANS
  3664.   if (eisnan (a)  || eisnan (b))
  3665.       return (-2);
  3666. #endif
  3667.   emovi (a, ai);
  3668.   p = ai;
  3669.   emovi (b, bi);
  3670.   q = bi;
  3671.  
  3672.   if (*p != *q)
  3673.     {                /* the signs are different */
  3674.       /* -0 equals + 0 */
  3675.       for (i = 1; i < NI - 1; i++)
  3676.     {
  3677.       if (ai[i] != 0)
  3678.         goto nzro;
  3679.       if (bi[i] != 0)
  3680.         goto nzro;
  3681.     }
  3682.       return (0);
  3683.     nzro:
  3684.       if (*p == 0)
  3685.     return (1);
  3686.       else
  3687.     return (-1);
  3688.     }
  3689.   /* both are the same sign */
  3690.   if (*p == 0)
  3691.     msign = 1;
  3692.   else
  3693.     msign = -1;
  3694.   i = NI - 1;
  3695.   do
  3696.     {
  3697.       if (*p++ != *q++)
  3698.     {
  3699.       goto diff;
  3700.     }
  3701.     }
  3702.   while (--i > 0);
  3703.  
  3704.   return (0);            /* equality */
  3705.  
  3706.  
  3707.  
  3708.  diff:
  3709.  
  3710.   if (*(--p) > *(--q))
  3711.     return (msign);        /* p is bigger */
  3712.   else
  3713.     return (-msign);        /* p is littler */
  3714. }
  3715.  
  3716.  
  3717.  
  3718.  
  3719. /* Find nearest integer to x = floor (x + 0.5)
  3720.  *
  3721.  * unsigned EMUSHORT x[NE], y[NE]
  3722.  * eround (x, y);
  3723.  */
  3724. void 
  3725. eround (x, y)
  3726.      unsigned EMUSHORT *x, *y;
  3727. {
  3728.   eadd (ehalf, x, y);
  3729.   efloor (y, y);
  3730. }
  3731.  
  3732.  
  3733.  
  3734.  
  3735. /*
  3736. ; convert HOST_WIDE_INT to e type
  3737. ;
  3738. ;    HOST_WIDE_INT l;
  3739. ;    unsigned EMUSHORT x[NE];
  3740. ;    ltoe (&l, x);
  3741. ; note &l is the memory address of l
  3742. */
  3743. void 
  3744. ltoe (lp, y)
  3745.      HOST_WIDE_INT *lp;
  3746.      unsigned EMUSHORT *y;
  3747. {
  3748.   unsigned EMUSHORT yi[NI];
  3749.   unsigned HOST_WIDE_INT ll;
  3750.   int k;
  3751.  
  3752.   ecleaz (yi);
  3753.   if (*lp < 0)
  3754.     {
  3755.       /* make it positive */
  3756.       ll = (unsigned HOST_WIDE_INT) (-(*lp));
  3757.       yi[0] = 0xffff;        /* put correct sign in the e type number */
  3758.     }
  3759.   else
  3760.     {
  3761.       ll = (unsigned HOST_WIDE_INT) (*lp);
  3762.     }
  3763.   /* move the long integer to yi significand area */
  3764. #if HOST_BITS_PER_WIDE_INT == 64
  3765.   yi[M] = (unsigned EMUSHORT) (ll >> 48);
  3766.   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
  3767.   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
  3768.   yi[M + 3] = (unsigned EMUSHORT) ll;
  3769.   yi[E] = EXONE + 47;        /* exponent if normalize shift count were 0 */
  3770. #else
  3771.   yi[M] = (unsigned EMUSHORT) (ll >> 16);
  3772.   yi[M + 1] = (unsigned EMUSHORT) ll;
  3773.   yi[E] = EXONE + 15;        /* exponent if normalize shift count were 0 */
  3774. #endif
  3775.  
  3776.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  3777.     ecleaz (yi);        /* it was zero */
  3778.   else
  3779.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  3780.   emovo (yi, y);        /* output the answer */
  3781. }
  3782.  
  3783. /*
  3784. ; convert unsigned HOST_WIDE_INT to e type
  3785. ;
  3786. ;    unsigned HOST_WIDE_INT l;
  3787. ;    unsigned EMUSHORT x[NE];
  3788. ;    ltox (&l, x);
  3789. ; note &l is the memory address of l
  3790. */
  3791. void 
  3792. ultoe (lp, y)
  3793.      unsigned HOST_WIDE_INT *lp;
  3794.      unsigned EMUSHORT *y;
  3795. {
  3796.   unsigned EMUSHORT yi[NI];
  3797.   unsigned HOST_WIDE_INT ll;
  3798.   int k;
  3799.  
  3800.   ecleaz (yi);
  3801.   ll = *lp;
  3802.  
  3803.   /* move the long integer to ayi significand area */
  3804. #if HOST_BITS_PER_WIDE_INT == 64
  3805.   yi[M] = (unsigned EMUSHORT) (ll >> 48);
  3806.   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
  3807.   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
  3808.   yi[M + 3] = (unsigned EMUSHORT) ll;
  3809.   yi[E] = EXONE + 47;        /* exponent if normalize shift count were 0 */
  3810. #else
  3811.   yi[M] = (unsigned EMUSHORT) (ll >> 16);
  3812.   yi[M + 1] = (unsigned EMUSHORT) ll;
  3813.   yi[E] = EXONE + 15;        /* exponent if normalize shift count were 0 */
  3814. #endif
  3815.  
  3816.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  3817.     ecleaz (yi);        /* it was zero */
  3818.   else
  3819.     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
  3820.   emovo (yi, y);        /* output the answer */
  3821. }
  3822.  
  3823.  
  3824. /*
  3825. ;    Find signed HOST_WIDE_INT integer and floating point fractional parts
  3826.  
  3827. ;    HOST_WIDE_INT i;
  3828. ;    unsigned EMUSHORT x[NE], frac[NE];
  3829. ;    xifrac (x, &i, frac);
  3830.  
  3831.   The integer output has the sign of the input.  The fraction is
  3832. the positive fractional part of abs (x).
  3833. */
  3834. void 
  3835. eifrac (x, i, frac)
  3836.      unsigned EMUSHORT *x;
  3837.      HOST_WIDE_INT *i;
  3838.      unsigned EMUSHORT *frac;
  3839. {
  3840.   unsigned EMUSHORT xi[NI];
  3841.   int j, k;
  3842.   unsigned HOST_WIDE_INT ll;
  3843.  
  3844.   emovi (x, xi);
  3845.   k = (int) xi[E] - (EXONE - 1);
  3846.   if (k <= 0)
  3847.     {
  3848.       /* if exponent <= 0, integer = 0 and real output is fraction */
  3849.       *i = 0L;
  3850.       emovo (xi, frac);
  3851.       return;
  3852.     }
  3853.   if (k > (HOST_BITS_PER_WIDE_INT - 1))
  3854.     {
  3855.       /* long integer overflow: output large integer
  3856.      and correct fraction  */
  3857.       if (xi[0])
  3858.     *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
  3859.       else
  3860.     *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
  3861.       eshift (xi, k);
  3862.       if (extra_warnings)
  3863.     warning ("overflow on truncation to integer");
  3864.     }
  3865.   else if (k > 16)
  3866.     {
  3867.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  3868.      then shift up by 16's.  */
  3869.       j = k - ((k >> 4) << 4);
  3870.       eshift (xi, j);
  3871.       ll = xi[M];
  3872.       k -= j;
  3873.       do
  3874.     {
  3875.       eshup6 (xi);
  3876.       ll = (ll << 16) | xi[M];
  3877.     }
  3878.       while ((k -= 16) > 0);
  3879.       *i = ll;
  3880.       if (xi[0])
  3881.     *i = -(*i);
  3882.     }
  3883.   else
  3884.       {
  3885.         /* shift not more than 16 bits */
  3886.           eshift (xi, k);
  3887.         *i = (HOST_WIDE_INT) xi[M] & 0xffff;
  3888.         if (xi[0])
  3889.       *i = -(*i);
  3890.       }
  3891.   xi[0] = 0;
  3892.   xi[E] = EXONE - 1;
  3893.   xi[M] = 0;
  3894.   if ((k = enormlz (xi)) > NBITS)
  3895.     ecleaz (xi);
  3896.   else
  3897.     xi[E] -= (unsigned EMUSHORT) k;
  3898.  
  3899.   emovo (xi, frac);
  3900. }
  3901.  
  3902.  
  3903. /* Find unsigned HOST_WIDE_INT integer and floating point fractional parts.
  3904.    A negative e type input yields integer output = 0
  3905.    but correct fraction.  */
  3906.  
  3907. void 
  3908. euifrac (x, i, frac)
  3909.      unsigned EMUSHORT *x;
  3910.      unsigned HOST_WIDE_INT *i;
  3911.      unsigned EMUSHORT *frac;
  3912. {
  3913.   unsigned HOST_WIDE_INT ll;
  3914.   unsigned EMUSHORT xi[NI];
  3915.   int j, k;
  3916.  
  3917.   emovi (x, xi);
  3918.   k = (int) xi[E] - (EXONE - 1);
  3919.   if (k <= 0)
  3920.     {
  3921.       /* if exponent <= 0, integer = 0 and argument is fraction */
  3922.       *i = 0L;
  3923.       emovo (xi, frac);
  3924.       return;
  3925.     }
  3926.   if (k > HOST_BITS_PER_WIDE_INT)
  3927.     {
  3928.       /* Long integer overflow: output large integer
  3929.      and correct fraction.
  3930.      Note, the BSD microvax compiler says that ~(0UL)
  3931.      is a syntax error.  */
  3932.       *i = ~(0L);
  3933.       eshift (xi, k);
  3934.       if (extra_warnings)
  3935.     warning ("overflow on truncation to unsigned integer");
  3936.     }
  3937.   else if (k > 16)
  3938.     {
  3939.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  3940.      then shift up by 16's.  */
  3941.       j = k - ((k >> 4) << 4);
  3942.       eshift (xi, j);
  3943.       ll = xi[M];
  3944.       k -= j;
  3945.       do
  3946.     {
  3947.       eshup6 (xi);
  3948.       ll = (ll << 16) | xi[M];
  3949.     }
  3950.       while ((k -= 16) > 0);
  3951.       *i = ll;
  3952.     }
  3953.   else
  3954.     {
  3955.       /* shift not more than 16 bits */
  3956.       eshift (xi, k);
  3957.       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
  3958.     }
  3959.  
  3960.   if (xi[0])  /* A negative value yields unsigned integer 0. */
  3961.     *i = 0L;
  3962.  
  3963.   xi[0] = 0;
  3964.   xi[E] = EXONE - 1;
  3965.   xi[M] = 0;
  3966.   if ((k = enormlz (xi)) > NBITS)
  3967.     ecleaz (xi);
  3968.   else
  3969.     xi[E] -= (unsigned EMUSHORT) k;
  3970.  
  3971.   emovo (xi, frac);
  3972. }
  3973.  
  3974.  
  3975.  
  3976. /*
  3977. ;    Shift significand
  3978. ;
  3979. ;    Shifts significand area up or down by the number of bits
  3980. ;    given by the variable sc.
  3981. */
  3982. int 
  3983. eshift (x, sc)
  3984.      unsigned EMUSHORT *x;
  3985.      int sc;
  3986. {
  3987.   unsigned EMUSHORT lost;
  3988.   unsigned EMUSHORT *p;
  3989.  
  3990.   if (sc == 0)
  3991.     return (0);
  3992.  
  3993.   lost = 0;
  3994.   p = x + NI - 1;
  3995.  
  3996.   if (sc < 0)
  3997.     {
  3998.       sc = -sc;
  3999.       while (sc >= 16)
  4000.     {
  4001.       lost |= *p;        /* remember lost bits */
  4002.       eshdn6 (x);
  4003.       sc -= 16;
  4004.     }
  4005.  
  4006.       while (sc >= 8)
  4007.     {
  4008.       lost |= *p & 0xff;
  4009.       eshdn8 (x);
  4010.       sc -= 8;
  4011.     }
  4012.  
  4013.       while (sc > 0)
  4014.     {
  4015.       lost |= *p & 1;
  4016.       eshdn1 (x);
  4017.       sc -= 1;
  4018.     }
  4019.     }
  4020.   else
  4021.     {
  4022.       while (sc >= 16)
  4023.     {
  4024.       eshup6 (x);
  4025.       sc -= 16;
  4026.     }
  4027.  
  4028.       while (sc >= 8)
  4029.     {
  4030.       eshup8 (x);
  4031.       sc -= 8;
  4032.     }
  4033.  
  4034.       while (sc > 0)
  4035.     {
  4036.       eshup1 (x);
  4037.       sc -= 1;
  4038.     }
  4039.     }
  4040.   if (lost)
  4041.     lost = 1;
  4042.   return ((int) lost);
  4043. }
  4044.  
  4045.  
  4046.  
  4047. /*
  4048. ;    normalize
  4049. ;
  4050. ; Shift normalizes the significand area pointed to by argument
  4051. ; shift count (up = positive) is returned.
  4052. */
  4053. int 
  4054. enormlz (x)
  4055.      unsigned EMUSHORT x[];
  4056. {
  4057.   register unsigned EMUSHORT *p;
  4058.   int sc;
  4059.  
  4060.   sc = 0;
  4061.   p = &x[M];
  4062.   if (*p != 0)
  4063.     goto normdn;
  4064.   ++p;
  4065.   if (*p & 0x8000)
  4066.     return (0);            /* already normalized */
  4067.   while (*p == 0)
  4068.     {
  4069.       eshup6 (x);
  4070.       sc += 16;
  4071.       /* With guard word, there are NBITS+16 bits available.
  4072.        * return true if all are zero.
  4073.        */
  4074.       if (sc > NBITS)
  4075.     return (sc);
  4076.     }
  4077.   /* see if high byte is zero */
  4078.   while ((*p & 0xff00) == 0)
  4079.     {
  4080.       eshup8 (x);
  4081.       sc += 8;
  4082.     }
  4083.   /* now shift 1 bit at a time */
  4084.   while ((*p & 0x8000) == 0)
  4085.     {
  4086.       eshup1 (x);
  4087.       sc += 1;
  4088.       if (sc > NBITS)
  4089.     {
  4090.       mtherr ("enormlz", UNDERFLOW);
  4091.       return (sc);
  4092.     }
  4093.     }
  4094.   return (sc);
  4095.  
  4096.   /* Normalize by shifting down out of the high guard word
  4097.      of the significand */
  4098.  normdn:
  4099.  
  4100.   if (*p & 0xff00)
  4101.     {
  4102.       eshdn8 (x);
  4103.       sc -= 8;
  4104.     }
  4105.   while (*p != 0)
  4106.     {
  4107.       eshdn1 (x);
  4108.       sc -= 1;
  4109.  
  4110.       if (sc < -NBITS)
  4111.     {
  4112.       mtherr ("enormlz", OVERFLOW);
  4113.       return (sc);
  4114.     }
  4115.     }
  4116.   return (sc);
  4117. }
  4118.  
  4119.  
  4120.  
  4121.  
  4122. /* Convert e type number to decimal format ASCII string.
  4123.  * The constants are for 64 bit precision.
  4124.  */
  4125.  
  4126. #define NTEN 12
  4127. #define MAXP 4096
  4128.  
  4129. #if LONG_DOUBLE_TYPE_SIZE == 128
  4130. static unsigned EMUSHORT etens[NTEN + 1][NE] =
  4131. {
  4132.   {0x6576, 0x4a92, 0x804a, 0x153f,
  4133.    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  4134.   {0x6a32, 0xce52, 0x329a, 0x28ce,
  4135.    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  4136.   {0x526c, 0x50ce, 0xf18b, 0x3d28,
  4137.    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  4138.   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
  4139.    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  4140.   {0x851e, 0xeab7, 0x98fe, 0x901b,
  4141.    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  4142.   {0x0235, 0x0137, 0x36b1, 0x336c,
  4143.    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  4144.   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
  4145.    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  4146.   {0x0000, 0x0000, 0x0000, 0x0000,
  4147.    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  4148.   {0x0000, 0x0000, 0x0000, 0x0000,
  4149.    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  4150.   {0x0000, 0x0000, 0x0000, 0x0000,
  4151.    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  4152.   {0x0000, 0x0000, 0x0000, 0x0000,
  4153.    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  4154.   {0x0000, 0x0000, 0x0000, 0x0000,
  4155.    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  4156.   {0x0000, 0x0000, 0x0000, 0x0000,
  4157.    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  4158. };
  4159.  
  4160. static unsigned EMUSHORT emtens[NTEN + 1][NE] =
  4161. {
  4162.   {0x2030, 0xcffc, 0xa1c3, 0x8123,
  4163.    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  4164.   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
  4165.    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  4166.   {0xf53f, 0xf698, 0x6bd3, 0x0158,
  4167.    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  4168.   {0xe731, 0x04d4, 0xe3f2, 0xd332,
  4169.    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  4170.   {0xa23e, 0x5308, 0xfefb, 0x1155,
  4171.    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  4172.   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
  4173.    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  4174.   {0x2a20, 0x6224, 0x47b3, 0x98d7,
  4175.    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  4176.   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
  4177.    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  4178.   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
  4179.    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  4180.   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
  4181.    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  4182.   {0xc155, 0xa4a8, 0x404e, 0x6113,
  4183.    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  4184.   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
  4185.    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  4186.   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
  4187.    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  4188. };
  4189. #else
  4190. /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
  4191. static unsigned EMUSHORT etens[NTEN + 1][NE] =
  4192. {
  4193.   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  4194.   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  4195.   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  4196.   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  4197.   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  4198.   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  4199.   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  4200.   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  4201.   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  4202.   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  4203.   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  4204.   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  4205.   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  4206. };
  4207.  
  4208. static unsigned EMUSHORT emtens[NTEN + 1][NE] =
  4209. {
  4210.   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  4211.   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  4212.   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  4213.   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  4214.   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  4215.   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  4216.   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  4217.   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  4218.   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  4219.   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  4220.   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  4221.   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  4222.   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  4223. };
  4224. #endif
  4225.  
  4226. void 
  4227. e24toasc (x, string, ndigs)
  4228.      unsigned EMUSHORT x[];
  4229.      char *string;
  4230.      int ndigs;
  4231. {
  4232.   unsigned EMUSHORT w[NI];
  4233.  
  4234.   e24toe (x, w);
  4235.   etoasc (w, string, ndigs);
  4236. }
  4237.  
  4238.  
  4239. void 
  4240. e53toasc (x, string, ndigs)
  4241.      unsigned EMUSHORT x[];
  4242.      char *string;
  4243.      int ndigs;
  4244. {
  4245.   unsigned EMUSHORT w[NI];
  4246.  
  4247.   e53toe (x, w);
  4248.   etoasc (w, string, ndigs);
  4249. }
  4250.  
  4251.  
  4252. void 
  4253. e64toasc (x, string, ndigs)
  4254.      unsigned EMUSHORT x[];
  4255.      char *string;
  4256.      int ndigs;
  4257. {
  4258.   unsigned EMUSHORT w[NI];
  4259.  
  4260.   e64toe (x, w);
  4261.   etoasc (w, string, ndigs);
  4262. }
  4263.  
  4264. void 
  4265. e113toasc (x, string, ndigs)
  4266.      unsigned EMUSHORT x[];
  4267.      char *string;
  4268.      int ndigs;
  4269. {
  4270.   unsigned EMUSHORT w[NI];
  4271.  
  4272.   e113toe (x, w);
  4273.   etoasc (w, string, ndigs);
  4274. }
  4275.  
  4276.  
  4277. static char wstring[80];    /* working storage for ASCII output */
  4278.  
  4279. void 
  4280. etoasc (x, string, ndigs)
  4281.      unsigned EMUSHORT x[];
  4282.      char *string;
  4283.      int ndigs;
  4284. {
  4285.   EMUSHORT digit;
  4286.   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
  4287.   unsigned EMUSHORT *p, *r, *ten;
  4288.   unsigned EMUSHORT sign;
  4289.   int i, j, k, expon, rndsav;
  4290.   char *s, *ss;
  4291.   unsigned EMUSHORT m;
  4292.  
  4293.  
  4294.   rndsav = rndprc;
  4295.   ss = string;
  4296.   s = wstring;
  4297.   *ss = '\0';
  4298.   *s = '\0';
  4299. #ifdef NANS
  4300.   if (eisnan (x))
  4301.     {
  4302.       sprintf (wstring, " NaN ");
  4303.       goto bxit;
  4304.     }
  4305. #endif
  4306.   rndprc = NBITS;        /* set to full precision */
  4307.   emov (x, y);            /* retain external format */
  4308.   if (y[NE - 1] & 0x8000)
  4309.     {
  4310.       sign = 0xffff;
  4311.       y[NE - 1] &= 0x7fff;
  4312.     }
  4313.   else
  4314.     {
  4315.       sign = 0;
  4316.     }
  4317.   expon = 0;
  4318.   ten = &etens[NTEN][0];
  4319.   emov (eone, t);
  4320.   /* Test for zero exponent */
  4321.   if (y[NE - 1] == 0)
  4322.     {
  4323.       for (k = 0; k < NE - 1; k++)
  4324.     {
  4325.       if (y[k] != 0)
  4326.         goto tnzro;        /* denormalized number */
  4327.     }
  4328.       goto isone;        /* legal all zeros */
  4329.     }
  4330.  tnzro:
  4331.  
  4332.   /* Test for infinity. */
  4333.   if (y[NE - 1] == 0x7fff)
  4334.     {
  4335.       if (sign)
  4336.     sprintf (wstring, " -Infinity ");
  4337.       else
  4338.     sprintf (wstring, " Infinity ");
  4339.       goto bxit;
  4340.     }
  4341.  
  4342.   /* Test for exponent nonzero but significand denormalized.
  4343.    * This is an error condition.
  4344.    */
  4345.   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
  4346.     {
  4347.       mtherr ("etoasc", DOMAIN);
  4348.       sprintf (wstring, "NaN");
  4349.       goto bxit;
  4350.     }
  4351.  
  4352.   /* Compare to 1.0 */
  4353.   i = ecmp (eone, y);
  4354.   if (i == 0)
  4355.     goto isone;
  4356.  
  4357.   if (i == -2)
  4358.     abort ();
  4359.  
  4360.   if (i < 0)
  4361.     {                /* Number is greater than 1 */
  4362.       /* Convert significand to an integer and strip trailing decimal zeros. */
  4363.       emov (y, u);
  4364.       u[NE - 1] = EXONE + NBITS - 1;
  4365.  
  4366.       p = &etens[NTEN - 4][0];
  4367.       m = 16;
  4368.       do
  4369.     {
  4370.       ediv (p, u, t);
  4371.       efloor (t, w);
  4372.       for (j = 0; j < NE - 1; j++)
  4373.         {
  4374.           if (t[j] != w[j])
  4375.         goto noint;
  4376.         }
  4377.       emov (t, u);
  4378.       expon += (int) m;
  4379.     noint:
  4380.       p += NE;
  4381.       m >>= 1;
  4382.     }
  4383.       while (m != 0);
  4384.  
  4385.       /* Rescale from integer significand */
  4386.       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
  4387.       emov (u, y);
  4388.       /* Find power of 10 */
  4389.       emov (eone, t);
  4390.       m = MAXP;
  4391.       p = &etens[0][0];
  4392.       /* An unordered compare result shouldn't happen here. */
  4393.       while (ecmp (ten, u) <= 0)
  4394.     {
  4395.       if (ecmp (p, u) <= 0)
  4396.         {
  4397.           ediv (p, u, u);
  4398.           emul (p, t, t);
  4399.           expon += (int) m;
  4400.         }
  4401.       m >>= 1;
  4402.       if (m == 0)
  4403.         break;
  4404.       p += NE;
  4405.     }
  4406.     }
  4407.   else
  4408.     {                /* Number is less than 1.0 */
  4409.       /* Pad significand with trailing decimal zeros. */
  4410.       if (y[NE - 1] == 0)
  4411.     {
  4412.       while ((y[NE - 2] & 0x8000) == 0)
  4413.         {
  4414.           emul (ten, y, y);
  4415.           expon -= 1;
  4416.         }
  4417.     }
  4418.       else
  4419.     {
  4420.       emovi (y, w);
  4421.       for (i = 0; i < NDEC + 1; i++)
  4422.         {
  4423.           if ((w[NI - 1] & 0x7) != 0)
  4424.         break;
  4425.           /* multiply by 10 */
  4426.           emovz (w, u);
  4427.           eshdn1 (u);
  4428.           eshdn1 (u);
  4429.           eaddm (w, u);
  4430.           u[1] += 3;
  4431.           while (u[2] != 0)
  4432.         {
  4433.           eshdn1 (u);
  4434.           u[1] += 1;
  4435.         }
  4436.           if (u[NI - 1] != 0)
  4437.         break;
  4438.           if (eone[NE - 1] <= u[1])
  4439.         break;
  4440.           emovz (u, w);
  4441.           expon -= 1;
  4442.         }
  4443.       emovo (w, y);
  4444.     }
  4445.       k = -MAXP;
  4446.       p = &emtens[0][0];
  4447.       r = &etens[0][0];
  4448.       emov (y, w);
  4449.       emov (eone, t);
  4450.       while (ecmp (eone, w) > 0)
  4451.     {
  4452.       if (ecmp (p, w) >= 0)
  4453.         {
  4454.           emul (r, w, w);
  4455.           emul (r, t, t);
  4456.           expon += k;
  4457.         }
  4458.       k /= 2;
  4459.       if (k == 0)
  4460.         break;
  4461.       p += NE;
  4462.       r += NE;
  4463.     }
  4464.       ediv (t, eone, t);
  4465.     }
  4466.  isone:
  4467.   /* Find the first (leading) digit. */
  4468.   emovi (t, w);
  4469.   emovz (w, t);
  4470.   emovi (y, w);
  4471.   emovz (w, y);
  4472.   eiremain (t, y);
  4473.   digit = equot[NI - 1];
  4474.   while ((digit == 0) && (ecmp (y, ezero) != 0))
  4475.     {
  4476.       eshup1 (y);
  4477.       emovz (y, u);
  4478.       eshup1 (u);
  4479.       eshup1 (u);
  4480.       eaddm (u, y);
  4481.       eiremain (t, y);
  4482.       digit = equot[NI - 1];
  4483.       expon -= 1;
  4484.     }
  4485.   s = wstring;
  4486.   if (sign)
  4487.     *s++ = '-';
  4488.   else
  4489.     *s++ = ' ';
  4490.   /* Examine number of digits requested by caller. */
  4491.   if (ndigs < 0)
  4492.     ndigs = 0;
  4493.   if (ndigs > NDEC)
  4494.     ndigs = NDEC;
  4495.   if (digit == 10)
  4496.     {
  4497.       *s++ = '1';
  4498.       *s++ = '.';
  4499.       if (ndigs > 0)
  4500.     {
  4501.       *s++ = '0';
  4502.       ndigs -= 1;
  4503.     }
  4504.       expon += 1;
  4505.     }
  4506.   else
  4507.     {
  4508.       *s++ = (char)digit + '0';
  4509.       *s++ = '.';
  4510.     }
  4511.   /* Generate digits after the decimal point. */
  4512.   for (k = 0; k <= ndigs; k++)
  4513.     {
  4514.       /* multiply current number by 10, without normalizing */
  4515.       eshup1 (y);
  4516.       emovz (y, u);
  4517.       eshup1 (u);
  4518.       eshup1 (u);
  4519.       eaddm (u, y);
  4520.       eiremain (t, y);
  4521.       *s++ = (char) equot[NI - 1] + '0';
  4522.     }
  4523.   digit = equot[NI - 1];
  4524.   --s;
  4525.   ss = s;
  4526.   /* round off the ASCII string */
  4527.   if (digit > 4)
  4528.     {
  4529.       /* Test for critical rounding case in ASCII output. */
  4530.       if (digit == 5)
  4531.     {
  4532.       emovo (y, t);
  4533.       if (ecmp (t, ezero) != 0)
  4534.         goto roun;        /* round to nearest */
  4535.       if ((*(s - 1) & 1) == 0)
  4536.         goto doexp;        /* round to even */
  4537.     }
  4538.       /* Round up and propagate carry-outs */
  4539.     roun:
  4540.       --s;
  4541.       k = *s & 0x7f;
  4542.       /* Carry out to most significant digit? */
  4543.       if (k == '.')
  4544.     {
  4545.       --s;
  4546.       k = *s;
  4547.       k += 1;
  4548.       *s = (char) k;
  4549.       /* Most significant digit carries to 10? */
  4550.       if (k > '9')
  4551.         {
  4552.           expon += 1;
  4553.           *s = '1';
  4554.         }
  4555.       goto doexp;
  4556.     }
  4557.       /* Round up and carry out from less significant digits */
  4558.       k += 1;
  4559.       *s = (char) k;
  4560.       if (k > '9')
  4561.     {
  4562.       *s = '0';
  4563.       goto roun;
  4564.     }
  4565.     }
  4566.  doexp:
  4567.   /*
  4568.      if (expon >= 0)
  4569.      sprintf (ss, "e+%d", expon);
  4570.      else
  4571.      sprintf (ss, "e%d", expon);
  4572.      */
  4573.   sprintf (ss, "e%d", expon);
  4574.  bxit:
  4575.   rndprc = rndsav;
  4576.   /* copy out the working string */
  4577.   s = string;
  4578.   ss = wstring;
  4579.   while (*ss == ' ')        /* strip possible leading space */
  4580.     ++ss;
  4581.   while ((*s++ = *ss++) != '\0')
  4582.     ;
  4583. }
  4584.  
  4585.  
  4586.  
  4587.  
  4588. /*
  4589. ;                                ASCTOQ
  4590. ;        ASCTOQ.MAC        LATEST REV: 11 JAN 84
  4591. ;                    SLM, 3 JAN 78
  4592. ;
  4593. ;    Convert ASCII string to quadruple precision floating point
  4594. ;
  4595. ;        Numeric input is free field decimal number
  4596. ;        with max of 15 digits with or without
  4597. ;        decimal point entered as ASCII from teletype.
  4598. ;    Entering E after the number followed by a second
  4599. ;    number causes the second number to be interpreted
  4600. ;    as a power of 10 to be multiplied by the first number
  4601. ;    (i.e., "scientific" notation).
  4602. ;
  4603. ;    Usage:
  4604. ;        asctoq (string, q);
  4605. */
  4606.  
  4607. /* ASCII to single */
  4608. void 
  4609. asctoe24 (s, y)
  4610.      char *s;
  4611.      unsigned EMUSHORT *y;
  4612. {
  4613.   asctoeg (s, y, 24);
  4614. }
  4615.  
  4616.  
  4617. /* ASCII to double */
  4618. void 
  4619. asctoe53 (s, y)
  4620.      char *s;
  4621.      unsigned EMUSHORT *y;
  4622. {
  4623. #if defined(DEC) || defined(IBM)
  4624.   asctoeg (s, y, 56);
  4625. #else
  4626.   asctoeg (s, y, 53);
  4627. #endif
  4628. }
  4629.  
  4630.  
  4631. /* ASCII to long double */
  4632. void 
  4633. asctoe64 (s, y)
  4634.      char *s;
  4635.      unsigned EMUSHORT *y;
  4636. {
  4637.   asctoeg (s, y, 64);
  4638. }
  4639.  
  4640. /* ASCII to 128-bit long double */
  4641. void 
  4642. asctoe113 (s, y)
  4643.      char *s;
  4644.      unsigned EMUSHORT *y;
  4645. {
  4646.   asctoeg (s, y, 113);
  4647. }
  4648.  
  4649. /* ASCII to super double */
  4650. void 
  4651. asctoe (s, y)
  4652.      char *s;
  4653.      unsigned EMUSHORT *y;
  4654. {
  4655.   asctoeg (s, y, NBITS);
  4656. }
  4657.  
  4658.  
  4659. /* ASCII to e type, with specified rounding precision = oprec. */
  4660. void 
  4661. asctoeg (ss, y, oprec)
  4662.      char *ss;
  4663.      unsigned EMUSHORT *y;
  4664.      int oprec;
  4665. {
  4666.   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
  4667.   int esign, decflg, sgnflg, nexp, exp, prec, lost;
  4668.   int k, trail, c, rndsav;
  4669.   EMULONG lexp;
  4670.   unsigned EMUSHORT nsign, *p;
  4671.   char *sp, *s, *lstr;
  4672.  
  4673.   /* Copy the input string. */
  4674.   lstr = (char *) alloca (strlen (ss) + 1);
  4675.   s = ss;
  4676.   while (*s == ' ')        /* skip leading spaces */
  4677.     ++s;
  4678.   sp = lstr;
  4679.   while ((*sp++ = *s++) != '\0')
  4680.     ;
  4681.   s = lstr;
  4682.  
  4683.   rndsav = rndprc;
  4684.   rndprc = NBITS;        /* Set to full precision */
  4685.   lost = 0;
  4686.   nsign = 0;
  4687.   decflg = 0;
  4688.   sgnflg = 0;
  4689.   nexp = 0;
  4690.   exp = 0;
  4691.   prec = 0;
  4692.   ecleaz (yy);
  4693.   trail = 0;
  4694.  
  4695.  nxtcom:
  4696.   k = *s - '0';
  4697.   if ((k >= 0) && (k <= 9))
  4698.     {
  4699.       /* Ignore leading zeros */
  4700.       if ((prec == 0) && (decflg == 0) && (k == 0))
  4701.     goto donchr;
  4702.       /* Identify and strip trailing zeros after the decimal point. */
  4703.       if ((trail == 0) && (decflg != 0))
  4704.     {
  4705.       sp = s;
  4706.       while ((*sp >= '0') && (*sp <= '9'))
  4707.         ++sp;
  4708.       /* Check for syntax error */
  4709.       c = *sp & 0x7f;
  4710.       if ((c != 'e') && (c != 'E') && (c != '\0')
  4711.           && (c != '\n') && (c != '\r') && (c != ' ')
  4712.           && (c != ','))
  4713.         goto error;
  4714.       --sp;
  4715.       while (*sp == '0')
  4716.         *sp-- = 'z';
  4717.       trail = 1;
  4718.       if (*s == 'z')
  4719.         goto donchr;
  4720.     }
  4721.       /* If enough digits were given to more than fill up the yy register,
  4722.        * continuing until overflow into the high guard word yy[2]
  4723.        * guarantees that there will be a roundoff bit at the top
  4724.        * of the low guard word after normalization.
  4725.        */
  4726.       if (yy[2] == 0)
  4727.     {
  4728.       if (decflg)
  4729.         nexp += 1;        /* count digits after decimal point */
  4730.       eshup1 (yy);        /* multiply current number by 10 */
  4731.       emovz (yy, xt);
  4732.       eshup1 (xt);
  4733.       eshup1 (xt);
  4734.       eaddm (xt, yy);
  4735.       ecleaz (xt);
  4736.       xt[NI - 2] = (unsigned EMUSHORT) k;
  4737.       eaddm (xt, yy);
  4738.     }
  4739.       else
  4740.     {
  4741.       /* Mark any lost non-zero digit.  */
  4742.       lost |= k;
  4743.       /* Count lost digits before the decimal point.  */
  4744.       if (decflg == 0)
  4745.         nexp -= 1;
  4746.     }
  4747.       prec += 1;
  4748.       goto donchr;
  4749.     }
  4750.  
  4751.   switch (*s)
  4752.     {
  4753.     case 'z':
  4754.       break;
  4755.     case 'E':
  4756.     case 'e':
  4757.       goto expnt;
  4758.     case '.':            /* decimal point */
  4759.       if (decflg)
  4760.     goto error;
  4761.       ++decflg;
  4762.       break;
  4763.     case '-':
  4764.       nsign = 0xffff;
  4765.       if (sgnflg)
  4766.     goto error;
  4767.       ++sgnflg;
  4768.       break;
  4769.     case '+':
  4770.       if (sgnflg)
  4771.     goto error;
  4772.       ++sgnflg;
  4773.       break;
  4774.     case ',':
  4775.     case ' ':
  4776.     case '\0':
  4777.     case '\n':
  4778.     case '\r':
  4779.       goto daldone;
  4780.     case 'i':
  4781.     case 'I':
  4782.       goto infinite;
  4783.     default:
  4784.     error:
  4785. #ifdef NANS
  4786.       einan (yy);
  4787. #else
  4788.       mtherr ("asctoe", DOMAIN);
  4789.       eclear (yy);
  4790. #endif
  4791.       goto aexit;
  4792.     }
  4793.  donchr:
  4794.   ++s;
  4795.   goto nxtcom;
  4796.  
  4797.   /* Exponent interpretation */
  4798.  expnt:
  4799.  
  4800.   esign = 1;
  4801.   exp = 0;
  4802.   ++s;
  4803.   /* check for + or - */
  4804.   if (*s == '-')
  4805.     {
  4806.       esign = -1;
  4807.       ++s;
  4808.     }
  4809.   if (*s == '+')
  4810.     ++s;
  4811.   while ((*s >= '0') && (*s <= '9'))
  4812.     {
  4813.       exp *= 10;
  4814.       exp += *s++ - '0';
  4815.       if (exp > -(MINDECEXP))
  4816.     {
  4817.       if (esign < 0)
  4818.         goto zero;
  4819.       else
  4820.         goto infinite;
  4821.     }
  4822.     }
  4823.   if (esign < 0)
  4824.     exp = -exp;
  4825.   if (exp > MAXDECEXP)
  4826.     {
  4827.  infinite:
  4828.       ecleaz (yy);
  4829.       yy[E] = 0x7fff;        /* infinity */
  4830.       goto aexit;
  4831.     }
  4832.   if (exp < MINDECEXP)
  4833.     {
  4834.  zero:
  4835.       ecleaz (yy);
  4836.       goto aexit;
  4837.     }
  4838.  
  4839.  daldone:
  4840.   nexp = exp - nexp;
  4841.   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
  4842.   while ((nexp > 0) && (yy[2] == 0))
  4843.     {
  4844.       emovz (yy, xt);
  4845.       eshup1 (xt);
  4846.       eshup1 (xt);
  4847.       eaddm (yy, xt);
  4848.       eshup1 (xt);
  4849.       if (xt[2] != 0)
  4850.     break;
  4851.       nexp -= 1;
  4852.       emovz (xt, yy);
  4853.     }
  4854.   if ((k = enormlz (yy)) > NBITS)
  4855.     {
  4856.       ecleaz (yy);
  4857.       goto aexit;
  4858.     }
  4859.   lexp = (EXONE - 1 + NBITS) - k;
  4860.   emdnorm (yy, lost, 0, lexp, 64);
  4861.   /* convert to external format */
  4862.  
  4863.  
  4864.   /* Multiply by 10**nexp.  If precision is 64 bits,
  4865.    * the maximum relative error incurred in forming 10**n
  4866.    * for 0 <= n <= 324 is 8.2e-20, at 10**180.
  4867.    * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
  4868.    * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
  4869.    */
  4870.   lexp = yy[E];
  4871.   if (nexp == 0)
  4872.     {
  4873.       k = 0;
  4874.       goto expdon;
  4875.     }
  4876.   esign = 1;
  4877.   if (nexp < 0)
  4878.     {
  4879.       nexp = -nexp;
  4880.       esign = -1;
  4881.       if (nexp > 4096)
  4882.     {            /* Punt.  Can't handle this without 2 divides. */
  4883.       emovi (etens[0], tt);
  4884.       lexp -= tt[E];
  4885.       k = edivm (tt, yy);
  4886.       lexp += EXONE;
  4887.       nexp -= 4096;
  4888.     }
  4889.     }
  4890.   p = &etens[NTEN][0];
  4891.   emov (eone, xt);
  4892.   exp = 1;
  4893.   do
  4894.     {
  4895.       if (exp & nexp)
  4896.     emul (p, xt, xt);
  4897.       p -= NE;
  4898.       exp = exp + exp;
  4899.     }
  4900.   while (exp <= MAXP);
  4901.  
  4902.   emovi (xt, tt);
  4903.   if (esign < 0)
  4904.     {
  4905.       lexp -= tt[E];
  4906.       k = edivm (tt, yy);
  4907.       lexp += EXONE;
  4908.     }
  4909.   else
  4910.     {
  4911.       lexp += tt[E];
  4912.       k = emulm (tt, yy);
  4913.       lexp -= EXONE - 1;
  4914.     }
  4915.  
  4916.  expdon:
  4917.  
  4918.   /* Round and convert directly to the destination type */
  4919.   if (oprec == 53)
  4920.     lexp -= EXONE - 0x3ff;
  4921. #ifdef IBM
  4922.   else if (oprec == 24 || oprec == 56)
  4923.     lexp -= EXONE - (0x41 << 2);
  4924. #else
  4925.   else if (oprec == 24)
  4926.     lexp -= EXONE - 0177;
  4927. #endif
  4928. #ifdef DEC
  4929.   else if (oprec == 56)
  4930.     lexp -= EXONE - 0201;
  4931. #endif
  4932.   rndprc = oprec;
  4933.   emdnorm (yy, k, 0, lexp, 64);
  4934.  
  4935.  aexit:
  4936.  
  4937.   rndprc = rndsav;
  4938.   yy[0] = nsign;
  4939.   switch (oprec)
  4940.     {
  4941. #ifdef DEC
  4942.     case 56:
  4943.       todec (yy, y);        /* see etodec.c */
  4944.       break;
  4945. #endif
  4946. #ifdef IBM
  4947.     case 56:
  4948.       toibm (yy, y, DFmode);
  4949.       break;
  4950. #endif
  4951.     case 53:
  4952.       toe53 (yy, y);
  4953.       break;
  4954.     case 24:
  4955.       toe24 (yy, y);
  4956.       break;
  4957.     case 64:
  4958.       toe64 (yy, y);
  4959.       break;
  4960.     case 113:
  4961.       toe113 (yy, y);
  4962.       break;
  4963.     case NBITS:
  4964.       emovo (yy, y);
  4965.       break;
  4966.     }
  4967. }
  4968.  
  4969.  
  4970.  
  4971. /* y = largest integer not greater than x
  4972.  * (truncated toward minus infinity)
  4973.  *
  4974.  * unsigned EMUSHORT x[NE], y[NE]
  4975.  *
  4976.  * efloor (x, y);
  4977.  */
  4978. static unsigned EMUSHORT bmask[] =
  4979. {
  4980.   0xffff,
  4981.   0xfffe,
  4982.   0xfffc,
  4983.   0xfff8,
  4984.   0xfff0,
  4985.   0xffe0,
  4986.   0xffc0,
  4987.   0xff80,
  4988.   0xff00,
  4989.   0xfe00,
  4990.   0xfc00,
  4991.   0xf800,
  4992.   0xf000,
  4993.   0xe000,
  4994.   0xc000,
  4995.   0x8000,
  4996.   0x0000,
  4997. };
  4998.  
  4999. void 
  5000. efloor (x, y)
  5001.      unsigned EMUSHORT x[], y[];
  5002. {
  5003.   register unsigned EMUSHORT *p;
  5004.   int e, expon, i;
  5005.   unsigned EMUSHORT f[NE];
  5006.  
  5007.   emov (x, f);            /* leave in external format */
  5008.   expon = (int) f[NE - 1];
  5009.   e = (expon & 0x7fff) - (EXONE - 1);
  5010.   if (e <= 0)
  5011.     {
  5012.       eclear (y);
  5013.       goto isitneg;
  5014.     }
  5015.   /* number of bits to clear out */
  5016.   e = NBITS - e;
  5017.   emov (f, y);
  5018.   if (e <= 0)
  5019.     return;
  5020.  
  5021.   p = &y[0];
  5022.   while (e >= 16)
  5023.     {
  5024.       *p++ = 0;
  5025.       e -= 16;
  5026.     }
  5027.   /* clear the remaining bits */
  5028.   *p &= bmask[e];
  5029.   /* truncate negatives toward minus infinity */
  5030.  isitneg:
  5031.  
  5032.   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
  5033.     {
  5034.       for (i = 0; i < NE - 1; i++)
  5035.     {
  5036.       if (f[i] != y[i])
  5037.         {
  5038.           esub (eone, y, y);
  5039.           break;
  5040.         }
  5041.     }
  5042.     }
  5043. }
  5044.  
  5045.  
  5046. /* unsigned EMUSHORT x[], s[];
  5047.  * int *exp;
  5048.  *
  5049.  * efrexp (x, exp, s);
  5050.  *
  5051.  * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
  5052.  * For example, 1.1 = 0.55 * 2**1
  5053.  * Handles denormalized numbers properly using long integer exp.
  5054.  */
  5055. void 
  5056. efrexp (x, exp, s)
  5057.      unsigned EMUSHORT x[];
  5058.      int *exp;
  5059.      unsigned EMUSHORT s[];
  5060. {
  5061.   unsigned EMUSHORT xi[NI];
  5062.   EMULONG li;
  5063.  
  5064.   emovi (x, xi);
  5065.   li = (EMULONG) ((EMUSHORT) xi[1]);
  5066.  
  5067.   if (li == 0)
  5068.     {
  5069.       li -= enormlz (xi);
  5070.     }
  5071.   xi[1] = 0x3ffe;
  5072.   emovo (xi, s);
  5073.   *exp = (int) (li - 0x3ffe);
  5074. }
  5075.  
  5076.  
  5077.  
  5078. /* unsigned EMUSHORT x[], y[];
  5079.  * int pwr2;
  5080.  *
  5081.  * eldexp (x, pwr2, y);
  5082.  *
  5083.  * Returns y = x * 2**pwr2.
  5084.  */
  5085. void 
  5086. eldexp (x, pwr2, y)
  5087.      unsigned EMUSHORT x[];
  5088.      int pwr2;
  5089.      unsigned EMUSHORT y[];
  5090. {
  5091.   unsigned EMUSHORT xi[NI];
  5092.   EMULONG li;
  5093.   int i;
  5094.  
  5095.   emovi (x, xi);
  5096.   li = xi[1];
  5097.   li += pwr2;
  5098.   i = 0;
  5099.   emdnorm (xi, i, i, li, 64);
  5100.   emovo (xi, y);
  5101. }
  5102.  
  5103.  
  5104. /* c = remainder after dividing b by a
  5105.  * Least significant integer quotient bits left in equot[].
  5106.  */
  5107. void 
  5108. eremain (a, b, c)
  5109.      unsigned EMUSHORT a[], b[], c[];
  5110. {
  5111.   unsigned EMUSHORT den[NI], num[NI];
  5112.  
  5113. #ifdef NANS
  5114.   if (eisinf (b)
  5115.       || (ecmp (a, ezero) == 0)
  5116.       || eisnan (a)
  5117.       || eisnan (b))
  5118.     {
  5119.       enan (c);
  5120.       return;
  5121.     }
  5122. #endif
  5123.   if (ecmp (a, ezero) == 0)
  5124.     {
  5125.       mtherr ("eremain", SING);
  5126.       eclear (c);
  5127.       return;
  5128.     }
  5129.   emovi (a, den);
  5130.   emovi (b, num);
  5131.   eiremain (den, num);
  5132.   /* Sign of remainder = sign of quotient */
  5133.   if (a[0] == b[0])
  5134.     num[0] = 0;
  5135.   else
  5136.     num[0] = 0xffff;
  5137.   emovo (num, c);
  5138. }
  5139.  
  5140. void 
  5141. eiremain (den, num)
  5142.      unsigned EMUSHORT den[], num[];
  5143. {
  5144.   EMULONG ld, ln;
  5145.   unsigned EMUSHORT j;
  5146.  
  5147.   ld = den[E];
  5148.   ld -= enormlz (den);
  5149.   ln = num[E];
  5150.   ln -= enormlz (num);
  5151.   ecleaz (equot);
  5152.   while (ln >= ld)
  5153.     {
  5154.       if (ecmpm (den, num) <= 0)
  5155.     {
  5156.       esubm (den, num);
  5157.       j = 1;
  5158.     }
  5159.       else
  5160.     {
  5161.       j = 0;
  5162.     }
  5163.       eshup1 (equot);
  5164.       equot[NI - 1] |= j;
  5165.       eshup1 (num);
  5166.       ln -= 1;
  5167.     }
  5168.   emdnorm (num, 0, 0, ln, 0);
  5169. }
  5170.  
  5171. /*                            mtherr.c
  5172.  *
  5173.  *    Library common error handling routine
  5174.  *
  5175.  *
  5176.  *
  5177.  * SYNOPSIS:
  5178.  *
  5179.  * char *fctnam;
  5180.  * int code;
  5181.  * void mtherr ();
  5182.  *
  5183.  * mtherr (fctnam, code);
  5184.  *
  5185.  *
  5186.  *
  5187.  * DESCRIPTION:
  5188.  *
  5189.  * This routine may be called to report one of the following
  5190.  * error conditions (in the include file mconf.h).
  5191.  *
  5192.  *   Mnemonic        Value          Significance
  5193.  *
  5194.  *    DOMAIN            1       argument domain error
  5195.  *    SING              2       function singularity
  5196.  *    OVERFLOW          3       overflow range error
  5197.  *    UNDERFLOW         4       underflow range error
  5198.  *    TLOSS             5       total loss of precision
  5199.  *    PLOSS             6       partial loss of precision
  5200.  *    INVALID           7       NaN - producing operation
  5201.  *    EDOM             33       Unix domain error code
  5202.  *    ERANGE           34       Unix range error code
  5203.  *
  5204.  * The default version of the file prints the function name,
  5205.  * passed to it by the pointer fctnam, followed by the
  5206.  * error condition.  The display is directed to the standard
  5207.  * output device.  The routine then returns to the calling
  5208.  * program.  Users may wish to modify the program to abort by
  5209.  * calling exit under severe error conditions such as domain
  5210.  * errors.
  5211.  *
  5212.  * Since all error conditions pass control to this function,
  5213.  * the display may be easily changed, eliminated, or directed
  5214.  * to an error logging device.
  5215.  *
  5216.  * SEE ALSO:
  5217.  *
  5218.  * mconf.h
  5219.  *
  5220.  */
  5221.  
  5222. /*
  5223. Cephes Math Library Release 2.0:  April, 1987
  5224. Copyright 1984, 1987 by Stephen L. Moshier
  5225. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  5226. */
  5227.  
  5228. /* include "mconf.h" */
  5229.  
  5230. /* Notice: the order of appearance of the following
  5231.  * messages is bound to the error codes defined
  5232.  * in mconf.h.
  5233.  */
  5234. #define NMSGS 8
  5235. static char *ermsg[NMSGS] =
  5236. {
  5237.   "unknown",            /* error code 0 */
  5238.   "domain",            /* error code 1 */
  5239.   "singularity",        /* et seq.      */
  5240.   "overflow",
  5241.   "underflow",
  5242.   "total loss of precision",
  5243.   "partial loss of precision",
  5244.   "invalid operation"
  5245. };
  5246.  
  5247. int merror = 0;
  5248. extern int merror;
  5249.  
  5250. void 
  5251. mtherr (name, code)
  5252.      char *name;
  5253.      int code;
  5254. {
  5255.   char errstr[80];
  5256.  
  5257.   /* Display string passed by calling program,
  5258.    * which is supposed to be the name of the
  5259.    * function in which the error occurred.
  5260.    */
  5261.  
  5262.   /* Display error message defined
  5263.    * by the code argument.
  5264.    */
  5265.   if ((code <= 0) || (code >= NMSGS))
  5266.     code = 0;
  5267.   sprintf (errstr, " %s %s error", name, ermsg[code]);
  5268.   if (extra_warnings)
  5269.     warning (errstr);
  5270.   /* Set global error message word */
  5271.   merror = code + 1;
  5272.  
  5273.   /* Return to calling
  5274.    * program
  5275.    */
  5276. }
  5277.  
  5278. #ifdef DEC
  5279. /* Here is etodec.c .
  5280.  *
  5281.  */
  5282.  
  5283. /*
  5284. ;    convert DEC double precision to e type
  5285. ;    double d;
  5286. ;    EMUSHORT e[NE];
  5287. ;    dectoe (&d, e);
  5288. */
  5289. void 
  5290. dectoe (d, e)
  5291.      unsigned EMUSHORT *d;
  5292.      unsigned EMUSHORT *e;
  5293. {
  5294.   unsigned EMUSHORT y[NI];
  5295.   register unsigned EMUSHORT r, *p;
  5296.  
  5297.   ecleaz (y);            /* start with a zero */
  5298.   p = y;            /* point to our number */
  5299.   r = *d;            /* get DEC exponent word */
  5300.   if (*d & (unsigned int) 0x8000)
  5301.     *p = 0xffff;        /* fill in our sign */
  5302.   ++p;                /* bump pointer to our exponent word */
  5303.   r &= 0x7fff;            /* strip the sign bit */
  5304.   if (r == 0)            /* answer = 0 if high order DEC word = 0 */
  5305.     goto done;
  5306.  
  5307.  
  5308.   r >>= 7;            /* shift exponent word down 7 bits */
  5309.   r += EXONE - 0201;        /* subtract DEC exponent offset */
  5310.   /* add our e type exponent offset */
  5311.   *p++ = r;            /* to form our exponent */
  5312.  
  5313.   r = *d++;            /* now do the high order mantissa */
  5314.   r &= 0177;            /* strip off the DEC exponent and sign bits */
  5315.   r |= 0200;            /* the DEC understood high order mantissa bit */
  5316.   *p++ = r;            /* put result in our high guard word */
  5317.  
  5318.   *p++ = *d++;            /* fill in the rest of our mantissa */
  5319.   *p++ = *d++;
  5320.   *p = *d;
  5321.  
  5322.   eshdn8 (y);            /* shift our mantissa down 8 bits */
  5323.  done:
  5324.   emovo (y, e);
  5325. }
  5326.  
  5327.  
  5328.  
  5329. /*
  5330. ;    convert e type to DEC double precision
  5331. ;    double d;
  5332. ;    EMUSHORT e[NE];
  5333. ;    etodec (e, &d);
  5334. */
  5335.  
  5336. void 
  5337. etodec (x, d)
  5338.      unsigned EMUSHORT *x, *d;
  5339. {
  5340.   unsigned EMUSHORT xi[NI];
  5341.   EMULONG exp;
  5342.   int rndsav;
  5343.  
  5344.   emovi (x, xi);
  5345.   exp = (EMULONG) xi[E] - (EXONE - 0201);    /* adjust exponent for offsets */
  5346. /* round off to nearest or even */
  5347.   rndsav = rndprc;
  5348.   rndprc = 56;
  5349.   emdnorm (xi, 0, 0, exp, 64);
  5350.   rndprc = rndsav;
  5351.   todec (xi, d);
  5352. }
  5353.  
  5354. void 
  5355. todec (x, y)
  5356.      unsigned EMUSHORT *x, *y;
  5357. {
  5358.   unsigned EMUSHORT i;
  5359.   unsigned EMUSHORT *p;
  5360.  
  5361.   p = x;
  5362.   *y = 0;
  5363.   if (*p++)
  5364.     *y = 0100000;
  5365.   i = *p++;
  5366.   if (i == 0)
  5367.     {
  5368.       *y++ = 0;
  5369.       *y++ = 0;
  5370.       *y++ = 0;
  5371.       *y++ = 0;
  5372.       return;
  5373.     }
  5374.   if (i > 0377)
  5375.     {
  5376.       *y++ |= 077777;
  5377.       *y++ = 0xffff;
  5378.       *y++ = 0xffff;
  5379.       *y++ = 0xffff;
  5380. #ifdef ERANGE
  5381.       errno = ERANGE;
  5382. #endif
  5383.       return;
  5384.     }
  5385.   i &= 0377;
  5386.   i <<= 7;
  5387.   eshup8 (x);
  5388.   x[M] &= 0177;
  5389.   i |= x[M];
  5390.   *y++ |= i;
  5391.   *y++ = x[M + 1];
  5392.   *y++ = x[M + 2];
  5393.   *y++ = x[M + 3];
  5394. }
  5395. #endif /* DEC */
  5396.  
  5397. #ifdef IBM
  5398. /* Here is etoibm
  5399.  *
  5400.  */
  5401.  
  5402. /*
  5403. ;    convert IBM single/double precision to e type
  5404. ;    single/double d;
  5405. ;    EMUSHORT e[NE];
  5406. ;    enum machine_mode mode;    SFmode/DFmode
  5407. ;    ibmtoe (&d, e, mode);
  5408. */
  5409. void 
  5410. ibmtoe (d, e, mode)
  5411.      unsigned EMUSHORT *d;
  5412.      unsigned EMUSHORT *e;
  5413.      enum machine_mode mode;
  5414. {
  5415.   unsigned EMUSHORT y[NI];
  5416.   register unsigned EMUSHORT r, *p;
  5417.   int rndsav;
  5418.  
  5419.   ecleaz (y);            /* start with a zero */
  5420.   p = y;            /* point to our number */
  5421.   r = *d;            /* get IBM exponent word */
  5422.   if (*d & (unsigned int) 0x8000)
  5423.     *p = 0xffff;        /* fill in our sign */
  5424.   ++p;                /* bump pointer to our exponent word */
  5425.   r &= 0x7f00;            /* strip the sign bit */
  5426.   r >>= 6;            /* shift exponent word down 6 bits */
  5427.                 /* in fact shift by 8 right and 2 left */
  5428.   r += EXONE - (0x41 << 2);    /* subtract IBM exponent offset */
  5429.                   /* add our e type exponent offset */
  5430.   *p++ = r;            /* to form our exponent */
  5431.  
  5432.   *p++ = *d++ & 0xff;        /* now do the high order mantissa */
  5433.                 /* strip off the IBM exponent and sign bits */
  5434.   if (mode != SFmode)        /* there are only 2 words in SFmode */
  5435.     {
  5436.       *p++ = *d++;        /* fill in the rest of our mantissa */
  5437.       *p++ = *d++;
  5438.     }
  5439.   *p = *d;
  5440.  
  5441.   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
  5442.     y[0] = y[E] = 0;
  5443.   else
  5444.     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
  5445.                   /* handle change in RADIX */
  5446.   emovo (y, e);
  5447. }
  5448.  
  5449.  
  5450.  
  5451. /*
  5452. ;    convert e type to IBM single/double precision
  5453. ;    single/double d;
  5454. ;    EMUSHORT e[NE];
  5455. ;    enum machine_mode mode;    SFmode/DFmode
  5456. ;    etoibm (e, &d, mode);
  5457. */
  5458.  
  5459. void 
  5460. etoibm (x, d, mode)
  5461.      unsigned EMUSHORT *x, *d;
  5462.      enum machine_mode mode;
  5463. {
  5464.   unsigned EMUSHORT xi[NI];
  5465.   EMULONG exp;
  5466.   int rndsav;
  5467.  
  5468.   emovi (x, xi);
  5469.   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));    /* adjust exponent for offsets */
  5470.                             /* round off to nearest or even */
  5471.   rndsav = rndprc;
  5472.   rndprc = 56;
  5473.   emdnorm (xi, 0, 0, exp, 64);
  5474.   rndprc = rndsav;
  5475.   toibm (xi, d, mode);
  5476. }
  5477.  
  5478. void 
  5479. toibm (x, y, mode)
  5480.      unsigned EMUSHORT *x, *y;
  5481.      enum machine_mode mode;
  5482. {
  5483.   unsigned EMUSHORT i;
  5484.   unsigned EMUSHORT *p;
  5485.   int r;
  5486.  
  5487.   p = x;
  5488.   *y = 0;
  5489.   if (*p++)
  5490.     *y = 0x8000;
  5491.   i = *p++;
  5492.   if (i == 0)
  5493.     {
  5494.       *y++ = 0;
  5495.       *y++ = 0;
  5496.       if (mode != SFmode)
  5497.     {
  5498.       *y++ = 0;
  5499.       *y++ = 0;
  5500.     }
  5501.       return;
  5502.     }
  5503.   r = i & 0x3;
  5504.   i >>= 2;
  5505.   if (i > 0x7f)
  5506.     {
  5507.       *y++ |= 0x7fff;
  5508.       *y++ = 0xffff;
  5509.       if (mode != SFmode)
  5510.     {
  5511.       *y++ = 0xffff;
  5512.       *y++ = 0xffff;
  5513.     }
  5514. #ifdef ERANGE
  5515.       errno = ERANGE;
  5516. #endif
  5517.       return;
  5518.     }
  5519.   i &= 0x7f;
  5520.   *y |= (i << 8);
  5521.   eshift (x, r + 5);
  5522.   *y++ |= x[M];
  5523.   *y++ = x[M + 1];
  5524.   if (mode != SFmode)
  5525.     {
  5526.       *y++ = x[M + 2];
  5527.       *y++ = x[M + 3];
  5528.     }
  5529. }
  5530. #endif /* IBM */
  5531.  
  5532. /* Output a binary NaN bit pattern in the target machine's format.  */
  5533.  
  5534. /* If special NaN bit patterns are required, define them in tm.h
  5535.    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
  5536.    patterns here. */
  5537. #ifdef TFMODE_NAN
  5538. TFMODE_NAN;
  5539. #else
  5540. #ifdef MIEEE
  5541. unsigned EMUSHORT TFnan[8] =
  5542.  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
  5543. #endif
  5544. #ifdef IBMPC
  5545. unsigned EMUSHORT TFnan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
  5546. #endif
  5547. #endif
  5548.  
  5549. #ifdef XFMODE_NAN
  5550. XFMODE_NAN;
  5551. #else
  5552. #ifdef MIEEE
  5553. unsigned EMUSHORT XFnan[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
  5554. #endif
  5555. #ifdef IBMPC
  5556. unsigned EMUSHORT XFnan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
  5557. #endif
  5558. #endif
  5559.  
  5560. #ifdef DFMODE_NAN
  5561. DFMODE_NAN;
  5562. #else
  5563. #ifdef MIEEE
  5564. unsigned EMUSHORT DFnan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
  5565. #endif
  5566. #ifdef IBMPC
  5567. unsigned EMUSHORT DFnan[4] = {0, 0, 0, 0xfff8};
  5568. #endif
  5569. #endif
  5570.  
  5571. #ifdef SFMODE_NAN
  5572. SFMODE_NAN;
  5573. #else
  5574. #ifdef MIEEE
  5575. unsigned EMUSHORT SFnan[2] = {0x7fff, 0xffff};
  5576. #endif
  5577. #ifdef IBMPC
  5578. unsigned EMUSHORT SFnan[2] = {0, 0xffc0};
  5579. #endif
  5580. #endif
  5581.  
  5582.  
  5583. void
  5584. make_nan (nan, mode)
  5585. unsigned EMUSHORT *nan;
  5586. enum machine_mode mode;
  5587. {
  5588.   int i, n;
  5589.   unsigned EMUSHORT *p;
  5590.  
  5591.   n = 0;
  5592.   switch (mode)
  5593.     {
  5594. /* Possibly the `reserved operand' patterns on a VAX can be
  5595.    used like NaN's, but probably not in the same way as IEEE. */
  5596. #if !defined(DEC) && !defined(IBM)
  5597.     case TFmode:
  5598.       n = 8;
  5599.       p = TFnan;
  5600.       break;
  5601.     case XFmode:
  5602.       n = 6;
  5603.       p = XFnan;
  5604.       break;
  5605.     case DFmode:
  5606.       n = 4;
  5607.       p = DFnan;
  5608.       break;
  5609.     case SFmode:
  5610.       n = 2;
  5611.       p = SFnan;
  5612.       break;
  5613. #endif
  5614.     default:
  5615.       abort ();
  5616.     }
  5617.   for (i=0; i < n; i++)
  5618.     *nan++ = *p++;
  5619. }
  5620.  
  5621. /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
  5622.    This is the inverse of the function `etarsingle' invoked by
  5623.    REAL_VALUE_TO_TARGET_SINGLE.  */
  5624.  
  5625. REAL_VALUE_TYPE
  5626. ereal_from_float (f)
  5627.      unsigned long f;
  5628. {
  5629.   REAL_VALUE_TYPE r;
  5630.   unsigned EMUSHORT s[2];
  5631.   unsigned EMUSHORT e[NE];
  5632.  
  5633.   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
  5634.    This is the inverse operation to what the function `endian' does.  */
  5635. #if FLOAT_WORDS_BIG_ENDIAN
  5636.   s[0] = (unsigned EMUSHORT) (f >> 16);
  5637.   s[1] = (unsigned EMUSHORT) f;
  5638. #else
  5639.   s[0] = (unsigned EMUSHORT) f;
  5640.   s[1] = (unsigned EMUSHORT) (f >> 16);
  5641. #endif
  5642.   /* Convert and promote the target float to E-type. */
  5643.   e24toe (s, e);
  5644.   /* Output E-type to REAL_VALUE_TYPE. */
  5645.   PUT_REAL (e, &r);
  5646.   return r;
  5647. }
  5648.  
  5649.  
  5650. /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
  5651.    This is the inverse of the function `etardouble' invoked by
  5652.    REAL_VALUE_TO_TARGET_DOUBLE.
  5653.  
  5654.    The DFmode is stored as an array of long ints
  5655.    with 32 bits of the value per each long.  The first element
  5656.    of the input array holds the bits that would come first in the
  5657.    target computer's memory.  */
  5658.  
  5659. REAL_VALUE_TYPE
  5660. ereal_from_double (d)
  5661.      unsigned long d[];
  5662. {
  5663.   REAL_VALUE_TYPE r;
  5664.   unsigned EMUSHORT s[4];
  5665.   unsigned EMUSHORT e[NE];
  5666.  
  5667.   /* Convert array of 32 bit pieces to equivalent array of 16 bit pieces.
  5668.    This is the inverse of `endian'.   */
  5669. #if FLOAT_WORDS_BIG_ENDIAN
  5670.   s[0] = (unsigned EMUSHORT) (d[0] >> 16);
  5671.   s[1] = (unsigned EMUSHORT) d[0];
  5672.   s[2] = (unsigned EMUSHORT) (d[1] >> 16);
  5673.   s[3] = (unsigned EMUSHORT) d[1];
  5674. #else
  5675.   s[0] = (unsigned EMUSHORT) d[0];
  5676.   s[1] = (unsigned EMUSHORT) (d[0] >> 16);
  5677.   s[2] = (unsigned EMUSHORT) d[1];
  5678.   s[3] = (unsigned EMUSHORT) (d[1] >> 16);
  5679. #endif
  5680.   /* Convert target double to E-type. */
  5681.   e53toe (s, e);
  5682.   /* Output E-type to REAL_VALUE_TYPE. */
  5683.   PUT_REAL (e, &r);
  5684.   return r;
  5685. }
  5686.  
  5687.  
  5688. /* Convert target computer unsigned 64-bit integer to e-type.
  5689.    The endian-ness of DImode follows the convention for integers,
  5690.    so we use WORDS_BIG_ENDIAN here, not FLOAT_WORDS_BIG_ENDIAN.  */
  5691.  
  5692. void
  5693. uditoe (di, e)
  5694.      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
  5695.      unsigned EMUSHORT *e;
  5696. {
  5697.   unsigned EMUSHORT yi[NI];
  5698.   int k;
  5699.  
  5700.   ecleaz (yi);
  5701. #if WORDS_BIG_ENDIAN
  5702.   for (k = M; k < M + 4; k++)
  5703.     yi[k] = *di++;
  5704. #else
  5705.   for (k = M + 3; k >= M; k--)
  5706.     yi[k] = *di++;
  5707. #endif
  5708.   yi[E] = EXONE + 47;    /* exponent if normalize shift count were 0 */
  5709.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  5710.     ecleaz (yi);        /* it was zero */
  5711.   else
  5712.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  5713.   emovo (yi, e);
  5714. }
  5715.  
  5716. /* Convert target computer signed 64-bit integer to e-type. */
  5717.  
  5718. void
  5719. ditoe (di, e)
  5720.      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
  5721.      unsigned EMUSHORT *e;
  5722. {
  5723.   unsigned EMULONG acc;
  5724.   unsigned EMUSHORT yi[NI];
  5725.   unsigned EMUSHORT carry;
  5726.   int k, sign;
  5727.  
  5728.   ecleaz (yi);
  5729. #if WORDS_BIG_ENDIAN
  5730.   for (k = M; k < M + 4; k++)
  5731.     yi[k] = *di++;
  5732. #else
  5733.   for (k = M + 3; k >= M; k--)
  5734.     yi[k] = *di++;
  5735. #endif
  5736.   /* Take absolute value */
  5737.   sign = 0;
  5738.   if (yi[M] & 0x8000)
  5739.     {
  5740.       sign = 1;
  5741.       carry = 0;
  5742.       for (k = M + 3; k >= M; k--)
  5743.     {
  5744.       acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
  5745.       yi[k] = acc;
  5746.       carry = 0;
  5747.       if (acc & 0x10000)
  5748.         carry = 1;
  5749.     }
  5750.     }
  5751.   yi[E] = EXONE + 47;    /* exponent if normalize shift count were 0 */
  5752.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  5753.     ecleaz (yi);        /* it was zero */
  5754.   else
  5755.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  5756.   emovo (yi, e);
  5757.   if (sign)
  5758.     eneg (e);
  5759. }
  5760.  
  5761.  
  5762. /* Convert e-type to unsigned 64-bit int. */
  5763.  
  5764. void 
  5765. etoudi (x, i)
  5766.      unsigned EMUSHORT *x;
  5767.      unsigned EMUSHORT *i;
  5768. {
  5769.   unsigned EMUSHORT xi[NI];
  5770.   int j, k;
  5771.  
  5772.   emovi (x, xi);
  5773.   if (xi[0])
  5774.     {
  5775.       xi[M] = 0;
  5776.       goto noshift;
  5777.     }
  5778.   k = (int) xi[E] - (EXONE - 1);
  5779.   if (k <= 0)
  5780.     {
  5781.       for (j = 0; j < 4; j++)
  5782.     *i++ = 0;
  5783.       return;
  5784.     }
  5785.   if (k > 64)
  5786.     {
  5787.       for (j = 0; j < 4; j++)
  5788.     *i++ = 0xffff;
  5789.       if (extra_warnings)
  5790.     warning ("overflow on truncation to integer");
  5791.       return;
  5792.     }
  5793.   if (k > 16)
  5794.     {
  5795.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  5796.      then shift up by 16's.  */
  5797.       j = k - ((k >> 4) << 4);
  5798.       if (j == 0)
  5799.     j = 16;
  5800.       eshift (xi, j);
  5801. #if WORDS_BIG_ENDIAN
  5802.       *i++ = xi[M];
  5803. #else
  5804.       i += 3;
  5805.       *i-- = xi[M];
  5806. #endif
  5807.       k -= j;
  5808.       do
  5809.     {
  5810.       eshup6 (xi);
  5811. #if WORDS_BIG_ENDIAN
  5812.       *i++ = xi[M];
  5813. #else
  5814.       *i-- = xi[M];
  5815. #endif
  5816.     }
  5817.       while ((k -= 16) > 0);
  5818.     }
  5819.   else
  5820.     {
  5821.         /* shift not more than 16 bits */
  5822.       eshift (xi, k);
  5823.  
  5824. noshift:
  5825.  
  5826. #if WORDS_BIG_ENDIAN
  5827.       i += 3;
  5828.       *i-- = xi[M];
  5829.       *i-- = 0;
  5830.       *i-- = 0;
  5831.       *i = 0;
  5832. #else
  5833.       *i++ = xi[M];
  5834.       *i++ = 0;
  5835.       *i++ = 0;
  5836.       *i = 0;
  5837. #endif
  5838.     }
  5839. }
  5840.  
  5841.  
  5842. /* Convert e-type to signed 64-bit int. */
  5843.  
  5844. void 
  5845. etodi (x, i)
  5846.      unsigned EMUSHORT *x;
  5847.      unsigned EMUSHORT *i;
  5848. {
  5849.   unsigned EMULONG acc;
  5850.   unsigned EMUSHORT xi[NI];
  5851.   unsigned EMUSHORT carry;
  5852.   unsigned EMUSHORT *isave;
  5853.   int j, k;
  5854.  
  5855.   emovi (x, xi);
  5856.   k = (int) xi[E] - (EXONE - 1);
  5857.   if (k <= 0)
  5858.     {
  5859.       for (j = 0; j < 4; j++)
  5860.     *i++ = 0;
  5861.       return;
  5862.     }
  5863.   if (k > 64)
  5864.     {
  5865.       for (j = 0; j < 4; j++)
  5866.     *i++ = 0xffff;
  5867.       if (extra_warnings)
  5868.     warning ("overflow on truncation to integer");
  5869.       return;
  5870.     }
  5871.   isave = i;
  5872.   if (k > 16)
  5873.     {
  5874.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  5875.      then shift up by 16's.  */
  5876.       j = k - ((k >> 4) << 4);
  5877.       if (j == 0)
  5878.     j = 16;
  5879.       eshift (xi, j);
  5880. #if WORDS_BIG_ENDIAN
  5881.       *i++ = xi[M];
  5882. #else
  5883.       i += 3;
  5884.       *i-- = xi[M];
  5885. #endif
  5886.       k -= j;
  5887.       do
  5888.     {
  5889.       eshup6 (xi);
  5890. #if WORDS_BIG_ENDIAN
  5891.       *i++ = xi[M];
  5892. #else
  5893.       *i-- = xi[M];
  5894. #endif
  5895.     }
  5896.       while ((k -= 16) > 0);
  5897.     }
  5898.   else
  5899.     {
  5900.         /* shift not more than 16 bits */
  5901.       eshift (xi, k);
  5902.  
  5903. #if WORDS_BIG_ENDIAN
  5904.       i += 3;
  5905.       *i = xi[M];
  5906.       *i-- = 0;
  5907.       *i-- = 0;
  5908.       *i = 0;
  5909. #else
  5910.       *i++ = xi[M];
  5911.       *i++ = 0;
  5912.       *i++ = 0;
  5913.       *i = 0;
  5914. #endif
  5915.     }
  5916.   /* Negate if negative */
  5917.   if (xi[0])
  5918.     {
  5919.       carry = 0;
  5920. #if WORDS_BIG_ENDIAN
  5921.       isave += 3;
  5922. #endif
  5923.       for (k = 0; k < 4; k++)
  5924.     {
  5925.       acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
  5926. #if WORDS_BIG_ENDIAN
  5927.       *isave-- = acc;
  5928. #else
  5929.       *isave++ = acc;
  5930. #endif
  5931.       carry = 0;
  5932.       if (acc & 0x10000)
  5933.         carry = 1;
  5934.     }
  5935.     }
  5936. }
  5937.  
  5938.  
  5939. /* Longhand square root routine. */
  5940.  
  5941.  
  5942. static int esqinited = 0;
  5943. static unsigned short sqrndbit[NI];
  5944.  
  5945. void 
  5946. esqrt (x, y)
  5947.      unsigned EMUSHORT *x, *y;
  5948. {
  5949.   unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
  5950.   EMULONG m, exp;
  5951.   int i, j, k, n, nlups;
  5952.  
  5953.   if (esqinited == 0)
  5954.     {
  5955.       ecleaz (sqrndbit);
  5956.       sqrndbit[NI - 2] = 1;
  5957.       esqinited = 1;
  5958.     }
  5959.   /* Check for arg <= 0 */
  5960.   i = ecmp (x, ezero);
  5961.   if (i <= 0)
  5962.     {
  5963. #ifdef NANS
  5964.       if (i == -2)
  5965.     {
  5966.       enan (y);
  5967.       return;
  5968.     }
  5969. #endif
  5970.       eclear (y);
  5971.       if (i < 0)
  5972.     mtherr ("esqrt", DOMAIN);
  5973.       return;
  5974.     }
  5975.  
  5976. #ifdef INFINITY
  5977.   if (eisinf (x))
  5978.     {
  5979.       eclear (y);
  5980.       einfin (y);
  5981.       return;
  5982.     }
  5983. #endif
  5984.   /* Bring in the arg and renormalize if it is denormal. */
  5985.   emovi (x, xx);
  5986.   m = (EMULONG) xx[1];        /* local long word exponent */
  5987.   if (m == 0)
  5988.     m -= enormlz (xx);
  5989.  
  5990.   /* Divide exponent by 2 */
  5991.   m -= 0x3ffe;
  5992.   exp = (unsigned short) ((m / 2) + 0x3ffe);
  5993.  
  5994.   /* Adjust if exponent odd */
  5995.   if ((m & 1) != 0)
  5996.     {
  5997.       if (m > 0)
  5998.     exp += 1;
  5999.       eshdn1 (xx);
  6000.     }
  6001.  
  6002.   ecleaz (sq);
  6003.   ecleaz (num);
  6004.   n = 8;            /* get 8 bits of result per inner loop */
  6005.   nlups = rndprc;
  6006.   j = 0;
  6007.  
  6008.   while (nlups > 0)
  6009.     {
  6010.       /* bring in next word of arg */
  6011.       if (j < NE)
  6012.     num[NI - 1] = xx[j + 3];
  6013.       /* Do additional bit on last outer loop, for roundoff. */
  6014.       if (nlups <= 8)
  6015.     n = nlups + 1;
  6016.       for (i = 0; i < n; i++)
  6017.     {
  6018.       /* Next 2 bits of arg */
  6019.       eshup1 (num);
  6020.       eshup1 (num);
  6021.       /* Shift up answer */
  6022.       eshup1 (sq);
  6023.       /* Make trial divisor */
  6024.       for (k = 0; k < NI; k++)
  6025.         temp[k] = sq[k];
  6026.       eshup1 (temp);
  6027.       eaddm (sqrndbit, temp);
  6028.       /* Subtract and insert answer bit if it goes in */
  6029.       if (ecmpm (temp, num) <= 0)
  6030.         {
  6031.           esubm (temp, num);
  6032.           sq[NI - 2] |= 1;
  6033.         }
  6034.     }
  6035.       nlups -= n;
  6036.       j += 1;
  6037.     }
  6038.  
  6039.   /* Adjust for extra, roundoff loop done. */
  6040.   exp += (NBITS - 1) - rndprc;
  6041.  
  6042.   /* Sticky bit = 1 if the remainder is nonzero. */
  6043.   k = 0;
  6044.   for (i = 3; i < NI; i++)
  6045.     k |= (int) num[i];
  6046.  
  6047.   /* Renormalize and round off. */
  6048.   emdnorm (sq, k, 0, exp, 64);
  6049.   emovo (sq, y);
  6050. }
  6051.  
  6052. #endif /* EMU_NON_COMPILE not defined */
  6053.